!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2015-2016  National Renewable Energy Laboratory
! Copyright (C) 2016-2018  Envision Energy USA, LTD
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!**********************************************************************************************************************************
!> This module contains the input/output parameters and routines for the BeamDyn module.
MODULE BeamDyn_IO

   USE BeamDyn_BldNdOuts_IO
   USE BeamDyn_Types
   USE BeamDyn_Subs
   USE NWTC_Library
   !USE NWTC_LAPACK

   IMPLICIT NONE

   TYPE(ProgDesc), PARAMETER:: BeamDyn_Ver = ProgDesc('BeamDyn', '','')




      ! Mass tolerance for checking the mass matrices:
   REAL(SiKi),    PARAMETER :: Mass0Tol = 1E-3                           ! Tolerance for checking mass values



! ===================================================================================================
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary.
! ===================================================================================================
! This code was generated by Write_ChckOutLst.m at 08-May-2025 15:51:30.


     ! Indices for computing output channels:
     ! NOTES:
     !    (1) These parameters are in the order stored in "OutListParameters.xlsx"
     !    (2) Array AllOuts() must be dimensioned to the value of the largest output parameter

     !  Time:

   INTEGER(IntKi), PARAMETER      :: Time           =   0


     ! Reaction Loads:

   INTEGER(IntKi), PARAMETER      :: RootFxr        =   1
   INTEGER(IntKi), PARAMETER      :: RootFyr        =   2
   INTEGER(IntKi), PARAMETER      :: RootFzr        =   3
   INTEGER(IntKi), PARAMETER      :: RootMxr        =   4
   INTEGER(IntKi), PARAMETER      :: RootMyr        =   5
   INTEGER(IntKi), PARAMETER      :: RootMzr        =   6


     ! Tip Motions:

   INTEGER(IntKi), PARAMETER      :: TipTDxr        =   7
   INTEGER(IntKi), PARAMETER      :: TipTDyr        =   8
   INTEGER(IntKi), PARAMETER      :: TipTDzr        =   9
   INTEGER(IntKi), PARAMETER      :: TipRDxr        =  10
   INTEGER(IntKi), PARAMETER      :: TipRDyr        =  11
   INTEGER(IntKi), PARAMETER      :: TipRDzr        =  12
   INTEGER(IntKi), PARAMETER      :: TipTVXg        =  13
   INTEGER(IntKi), PARAMETER      :: TipTVYg        =  14
   INTEGER(IntKi), PARAMETER      :: TipTVZg        =  15
   INTEGER(IntKi), PARAMETER      :: TipRVXg        =  16
   INTEGER(IntKi), PARAMETER      :: TipRVYg        =  17
   INTEGER(IntKi), PARAMETER      :: TipRVZg        =  18
   INTEGER(IntKi), PARAMETER      :: TipTAXl        =  19
   INTEGER(IntKi), PARAMETER      :: TipTAYl        =  20
   INTEGER(IntKi), PARAMETER      :: TipTAZl        =  21
   INTEGER(IntKi), PARAMETER      :: TipRAXl        =  22
   INTEGER(IntKi), PARAMETER      :: TipRAYl        =  23
   INTEGER(IntKi), PARAMETER      :: TipRAZl        =  24


     ! Sectional Loads:

   INTEGER(IntKi), PARAMETER      :: N1Fxl          =  25
   INTEGER(IntKi), PARAMETER      :: N1Fyl          =  26
   INTEGER(IntKi), PARAMETER      :: N1Fzl          =  27
   INTEGER(IntKi), PARAMETER      :: N2Fxl          =  28
   INTEGER(IntKi), PARAMETER      :: N2Fyl          =  29
   INTEGER(IntKi), PARAMETER      :: N2Fzl          =  30
   INTEGER(IntKi), PARAMETER      :: N3Fxl          =  31
   INTEGER(IntKi), PARAMETER      :: N3Fyl          =  32
   INTEGER(IntKi), PARAMETER      :: N3Fzl          =  33
   INTEGER(IntKi), PARAMETER      :: N4Fxl          =  34
   INTEGER(IntKi), PARAMETER      :: N4Fyl          =  35
   INTEGER(IntKi), PARAMETER      :: N4Fzl          =  36
   INTEGER(IntKi), PARAMETER      :: N5Fxl          =  37
   INTEGER(IntKi), PARAMETER      :: N5Fyl          =  38
   INTEGER(IntKi), PARAMETER      :: N5Fzl          =  39
   INTEGER(IntKi), PARAMETER      :: N6Fxl          =  40
   INTEGER(IntKi), PARAMETER      :: N6Fyl          =  41
   INTEGER(IntKi), PARAMETER      :: N6Fzl          =  42
   INTEGER(IntKi), PARAMETER      :: N7Fxl          =  43
   INTEGER(IntKi), PARAMETER      :: N7Fyl          =  44
   INTEGER(IntKi), PARAMETER      :: N7Fzl          =  45
   INTEGER(IntKi), PARAMETER      :: N8Fxl          =  46
   INTEGER(IntKi), PARAMETER      :: N8Fyl          =  47
   INTEGER(IntKi), PARAMETER      :: N8Fzl          =  48
   INTEGER(IntKi), PARAMETER      :: N9Fxl          =  49
   INTEGER(IntKi), PARAMETER      :: N9Fyl          =  50
   INTEGER(IntKi), PARAMETER      :: N9Fzl          =  51
   INTEGER(IntKi), PARAMETER      :: N1Mxl          =  52
   INTEGER(IntKi), PARAMETER      :: N1Myl          =  53
   INTEGER(IntKi), PARAMETER      :: N1Mzl          =  54
   INTEGER(IntKi), PARAMETER      :: N2Mxl          =  55
   INTEGER(IntKi), PARAMETER      :: N2Myl          =  56
   INTEGER(IntKi), PARAMETER      :: N2Mzl          =  57
   INTEGER(IntKi), PARAMETER      :: N3Mxl          =  58
   INTEGER(IntKi), PARAMETER      :: N3Myl          =  59
   INTEGER(IntKi), PARAMETER      :: N3Mzl          =  60
   INTEGER(IntKi), PARAMETER      :: N4Mxl          =  61
   INTEGER(IntKi), PARAMETER      :: N4Myl          =  62
   INTEGER(IntKi), PARAMETER      :: N4Mzl          =  63
   INTEGER(IntKi), PARAMETER      :: N5Mxl          =  64
   INTEGER(IntKi), PARAMETER      :: N5Myl          =  65
   INTEGER(IntKi), PARAMETER      :: N5Mzl          =  66
   INTEGER(IntKi), PARAMETER      :: N6Mxl          =  67
   INTEGER(IntKi), PARAMETER      :: N6Myl          =  68
   INTEGER(IntKi), PARAMETER      :: N6Mzl          =  69
   INTEGER(IntKi), PARAMETER      :: N7Mxl          =  70
   INTEGER(IntKi), PARAMETER      :: N7Myl          =  71
   INTEGER(IntKi), PARAMETER      :: N7Mzl          =  72
   INTEGER(IntKi), PARAMETER      :: N8Mxl          =  73
   INTEGER(IntKi), PARAMETER      :: N8Myl          =  74
   INTEGER(IntKi), PARAMETER      :: N8Mzl          =  75
   INTEGER(IntKi), PARAMETER      :: N9Mxl          =  76
   INTEGER(IntKi), PARAMETER      :: N9Myl          =  77
   INTEGER(IntKi), PARAMETER      :: N9Mzl          =  78


     ! Sectional Motions:

   INTEGER(IntKi), PARAMETER      :: N1TDxr         =  79
   INTEGER(IntKi), PARAMETER      :: N1TDyr         =  80
   INTEGER(IntKi), PARAMETER      :: N1TDzr         =  81
   INTEGER(IntKi), PARAMETER      :: N2TDxr         =  82
   INTEGER(IntKi), PARAMETER      :: N2TDyr         =  83
   INTEGER(IntKi), PARAMETER      :: N2TDzr         =  84
   INTEGER(IntKi), PARAMETER      :: N3TDxr         =  85
   INTEGER(IntKi), PARAMETER      :: N3TDyr         =  86
   INTEGER(IntKi), PARAMETER      :: N3TDzr         =  87
   INTEGER(IntKi), PARAMETER      :: N4TDxr         =  88
   INTEGER(IntKi), PARAMETER      :: N4TDyr         =  89
   INTEGER(IntKi), PARAMETER      :: N4TDzr         =  90
   INTEGER(IntKi), PARAMETER      :: N5TDxr         =  91
   INTEGER(IntKi), PARAMETER      :: N5TDyr         =  92
   INTEGER(IntKi), PARAMETER      :: N5TDzr         =  93
   INTEGER(IntKi), PARAMETER      :: N6TDxr         =  94
   INTEGER(IntKi), PARAMETER      :: N6TDyr         =  95
   INTEGER(IntKi), PARAMETER      :: N6TDzr         =  96
   INTEGER(IntKi), PARAMETER      :: N7TDxr         =  97
   INTEGER(IntKi), PARAMETER      :: N7TDyr         =  98
   INTEGER(IntKi), PARAMETER      :: N7TDzr         =  99
   INTEGER(IntKi), PARAMETER      :: N8TDxr         = 100
   INTEGER(IntKi), PARAMETER      :: N8TDyr         = 101
   INTEGER(IntKi), PARAMETER      :: N8TDzr         = 102
   INTEGER(IntKi), PARAMETER      :: N9TDxr         = 103
   INTEGER(IntKi), PARAMETER      :: N9TDyr         = 104
   INTEGER(IntKi), PARAMETER      :: N9TDzr         = 105
   INTEGER(IntKi), PARAMETER      :: N1RDxr         = 106
   INTEGER(IntKi), PARAMETER      :: N1RDyr         = 107
   INTEGER(IntKi), PARAMETER      :: N1RDzr         = 108
   INTEGER(IntKi), PARAMETER      :: N2RDxr         = 109
   INTEGER(IntKi), PARAMETER      :: N2RDyr         = 110
   INTEGER(IntKi), PARAMETER      :: N2RDzr         = 111
   INTEGER(IntKi), PARAMETER      :: N3RDxr         = 112
   INTEGER(IntKi), PARAMETER      :: N3RDyr         = 113
   INTEGER(IntKi), PARAMETER      :: N3RDzr         = 114
   INTEGER(IntKi), PARAMETER      :: N4RDxr         = 115
   INTEGER(IntKi), PARAMETER      :: N4RDyr         = 116
   INTEGER(IntKi), PARAMETER      :: N4RDzr         = 117
   INTEGER(IntKi), PARAMETER      :: N5RDxr         = 118
   INTEGER(IntKi), PARAMETER      :: N5RDyr         = 119
   INTEGER(IntKi), PARAMETER      :: N5RDzr         = 120
   INTEGER(IntKi), PARAMETER      :: N6RDxr         = 121
   INTEGER(IntKi), PARAMETER      :: N6RDyr         = 122
   INTEGER(IntKi), PARAMETER      :: N6RDzr         = 123
   INTEGER(IntKi), PARAMETER      :: N7RDxr         = 124
   INTEGER(IntKi), PARAMETER      :: N7RDyr         = 125
   INTEGER(IntKi), PARAMETER      :: N7RDzr         = 126
   INTEGER(IntKi), PARAMETER      :: N8RDxr         = 127
   INTEGER(IntKi), PARAMETER      :: N8RDyr         = 128
   INTEGER(IntKi), PARAMETER      :: N8RDzr         = 129
   INTEGER(IntKi), PARAMETER      :: N9RDxr         = 130
   INTEGER(IntKi), PARAMETER      :: N9RDyr         = 131
   INTEGER(IntKi), PARAMETER      :: N9RDzr         = 132
   INTEGER(IntKi), PARAMETER      :: N1TVXg         = 133
   INTEGER(IntKi), PARAMETER      :: N1TVYg         = 134
   INTEGER(IntKi), PARAMETER      :: N1TVZg         = 135
   INTEGER(IntKi), PARAMETER      :: N2TVXg         = 136
   INTEGER(IntKi), PARAMETER      :: N2TVYg         = 137
   INTEGER(IntKi), PARAMETER      :: N2TVZg         = 138
   INTEGER(IntKi), PARAMETER      :: N3TVXg         = 139
   INTEGER(IntKi), PARAMETER      :: N3TVYg         = 140
   INTEGER(IntKi), PARAMETER      :: N3TVZg         = 141
   INTEGER(IntKi), PARAMETER      :: N4TVXg         = 142
   INTEGER(IntKi), PARAMETER      :: N4TVYg         = 143
   INTEGER(IntKi), PARAMETER      :: N4TVZg         = 144
   INTEGER(IntKi), PARAMETER      :: N5TVXg         = 145
   INTEGER(IntKi), PARAMETER      :: N5TVYg         = 146
   INTEGER(IntKi), PARAMETER      :: N5TVZg         = 147
   INTEGER(IntKi), PARAMETER      :: N6TVXg         = 148
   INTEGER(IntKi), PARAMETER      :: N6TVYg         = 149
   INTEGER(IntKi), PARAMETER      :: N6TVZg         = 150
   INTEGER(IntKi), PARAMETER      :: N7TVXg         = 151
   INTEGER(IntKi), PARAMETER      :: N7TVYg         = 152
   INTEGER(IntKi), PARAMETER      :: N7TVZg         = 153
   INTEGER(IntKi), PARAMETER      :: N8TVXg         = 154
   INTEGER(IntKi), PARAMETER      :: N8TVYg         = 155
   INTEGER(IntKi), PARAMETER      :: N8TVZg         = 156
   INTEGER(IntKi), PARAMETER      :: N9TVXg         = 157
   INTEGER(IntKi), PARAMETER      :: N9TVYg         = 158
   INTEGER(IntKi), PARAMETER      :: N9TVZg         = 159
   INTEGER(IntKi), PARAMETER      :: N1RVXg         = 160
   INTEGER(IntKi), PARAMETER      :: N1RVYg         = 161
   INTEGER(IntKi), PARAMETER      :: N1RVZg         = 162
   INTEGER(IntKi), PARAMETER      :: N2RVXg         = 163
   INTEGER(IntKi), PARAMETER      :: N2RVYg         = 164
   INTEGER(IntKi), PARAMETER      :: N2RVZg         = 165
   INTEGER(IntKi), PARAMETER      :: N3RVXg         = 166
   INTEGER(IntKi), PARAMETER      :: N3RVYg         = 167
   INTEGER(IntKi), PARAMETER      :: N3RVZg         = 168
   INTEGER(IntKi), PARAMETER      :: N4RVXg         = 169
   INTEGER(IntKi), PARAMETER      :: N4RVYg         = 170
   INTEGER(IntKi), PARAMETER      :: N4RVZg         = 171
   INTEGER(IntKi), PARAMETER      :: N5RVXg         = 172
   INTEGER(IntKi), PARAMETER      :: N5RVYg         = 173
   INTEGER(IntKi), PARAMETER      :: N5RVZg         = 174
   INTEGER(IntKi), PARAMETER      :: N6RVXg         = 175
   INTEGER(IntKi), PARAMETER      :: N6RVYg         = 176
   INTEGER(IntKi), PARAMETER      :: N6RVZg         = 177
   INTEGER(IntKi), PARAMETER      :: N7RVXg         = 178
   INTEGER(IntKi), PARAMETER      :: N7RVYg         = 179
   INTEGER(IntKi), PARAMETER      :: N7RVZg         = 180
   INTEGER(IntKi), PARAMETER      :: N8RVXg         = 181
   INTEGER(IntKi), PARAMETER      :: N8RVYg         = 182
   INTEGER(IntKi), PARAMETER      :: N8RVZg         = 183
   INTEGER(IntKi), PARAMETER      :: N9RVXg         = 184
   INTEGER(IntKi), PARAMETER      :: N9RVYg         = 185
   INTEGER(IntKi), PARAMETER      :: N9RVZg         = 186
   INTEGER(IntKi), PARAMETER      :: N1TAXl         = 187
   INTEGER(IntKi), PARAMETER      :: N1TAYl         = 188
   INTEGER(IntKi), PARAMETER      :: N1TAZl         = 189
   INTEGER(IntKi), PARAMETER      :: N2TAXl         = 190
   INTEGER(IntKi), PARAMETER      :: N2TAYl         = 191
   INTEGER(IntKi), PARAMETER      :: N2TAZl         = 192
   INTEGER(IntKi), PARAMETER      :: N3TAXl         = 193
   INTEGER(IntKi), PARAMETER      :: N3TAYl         = 194
   INTEGER(IntKi), PARAMETER      :: N3TAZl         = 195
   INTEGER(IntKi), PARAMETER      :: N4TAXl         = 196
   INTEGER(IntKi), PARAMETER      :: N4TAYl         = 197
   INTEGER(IntKi), PARAMETER      :: N4TAZl         = 198
   INTEGER(IntKi), PARAMETER      :: N5TAXl         = 199
   INTEGER(IntKi), PARAMETER      :: N5TAYl         = 200
   INTEGER(IntKi), PARAMETER      :: N5TAZl         = 201
   INTEGER(IntKi), PARAMETER      :: N6TAXl         = 202
   INTEGER(IntKi), PARAMETER      :: N6TAYl         = 203
   INTEGER(IntKi), PARAMETER      :: N6TAZl         = 204
   INTEGER(IntKi), PARAMETER      :: N7TAXl         = 205
   INTEGER(IntKi), PARAMETER      :: N7TAYl         = 206
   INTEGER(IntKi), PARAMETER      :: N7TAZl         = 207
   INTEGER(IntKi), PARAMETER      :: N8TAXl         = 208
   INTEGER(IntKi), PARAMETER      :: N8TAYl         = 209
   INTEGER(IntKi), PARAMETER      :: N8TAZl         = 210
   INTEGER(IntKi), PARAMETER      :: N9TAXl         = 211
   INTEGER(IntKi), PARAMETER      :: N9TAYl         = 212
   INTEGER(IntKi), PARAMETER      :: N9TAZl         = 213
   INTEGER(IntKi), PARAMETER      :: N1RAXl         = 214
   INTEGER(IntKi), PARAMETER      :: N1RAYl         = 215
   INTEGER(IntKi), PARAMETER      :: N1RAZl         = 216
   INTEGER(IntKi), PARAMETER      :: N2RAXl         = 217
   INTEGER(IntKi), PARAMETER      :: N2RAYl         = 218
   INTEGER(IntKi), PARAMETER      :: N2RAZl         = 219
   INTEGER(IntKi), PARAMETER      :: N3RAXl         = 220
   INTEGER(IntKi), PARAMETER      :: N3RAYl         = 221
   INTEGER(IntKi), PARAMETER      :: N3RAZl         = 222
   INTEGER(IntKi), PARAMETER      :: N4RAXl         = 223
   INTEGER(IntKi), PARAMETER      :: N4RAYl         = 224
   INTEGER(IntKi), PARAMETER      :: N4RAZl         = 225
   INTEGER(IntKi), PARAMETER      :: N5RAXl         = 226
   INTEGER(IntKi), PARAMETER      :: N5RAYl         = 227
   INTEGER(IntKi), PARAMETER      :: N5RAZl         = 228
   INTEGER(IntKi), PARAMETER      :: N6RAXl         = 229
   INTEGER(IntKi), PARAMETER      :: N6RAYl         = 230
   INTEGER(IntKi), PARAMETER      :: N6RAZl         = 231
   INTEGER(IntKi), PARAMETER      :: N7RAXl         = 232
   INTEGER(IntKi), PARAMETER      :: N7RAYl         = 233
   INTEGER(IntKi), PARAMETER      :: N7RAZl         = 234
   INTEGER(IntKi), PARAMETER      :: N8RAXl         = 235
   INTEGER(IntKi), PARAMETER      :: N8RAYl         = 236
   INTEGER(IntKi), PARAMETER      :: N8RAZl         = 237
   INTEGER(IntKi), PARAMETER      :: N9RAXl         = 238
   INTEGER(IntKi), PARAMETER      :: N9RAYl         = 239
   INTEGER(IntKi), PARAMETER      :: N9RAZl         = 240


     ! Applied Loads:

   INTEGER(IntKi), PARAMETER      :: N1PFxl         = 241
   INTEGER(IntKi), PARAMETER      :: N1PFyl         = 242
   INTEGER(IntKi), PARAMETER      :: N1PFzl         = 243
   INTEGER(IntKi), PARAMETER      :: N2PFxl         = 244
   INTEGER(IntKi), PARAMETER      :: N2PFyl         = 245
   INTEGER(IntKi), PARAMETER      :: N2PFzl         = 246
   INTEGER(IntKi), PARAMETER      :: N3PFxl         = 247
   INTEGER(IntKi), PARAMETER      :: N3PFyl         = 248
   INTEGER(IntKi), PARAMETER      :: N3PFzl         = 249
   INTEGER(IntKi), PARAMETER      :: N4PFxl         = 250
   INTEGER(IntKi), PARAMETER      :: N4PFyl         = 251
   INTEGER(IntKi), PARAMETER      :: N4PFzl         = 252
   INTEGER(IntKi), PARAMETER      :: N5PFxl         = 253
   INTEGER(IntKi), PARAMETER      :: N5PFyl         = 254
   INTEGER(IntKi), PARAMETER      :: N5PFzl         = 255
   INTEGER(IntKi), PARAMETER      :: N6PFxl         = 256
   INTEGER(IntKi), PARAMETER      :: N6PFyl         = 257
   INTEGER(IntKi), PARAMETER      :: N6PFzl         = 258
   INTEGER(IntKi), PARAMETER      :: N7PFxl         = 259
   INTEGER(IntKi), PARAMETER      :: N7PFyl         = 260
   INTEGER(IntKi), PARAMETER      :: N7PFzl         = 261
   INTEGER(IntKi), PARAMETER      :: N8PFxl         = 262
   INTEGER(IntKi), PARAMETER      :: N8PFyl         = 263
   INTEGER(IntKi), PARAMETER      :: N8PFzl         = 264
   INTEGER(IntKi), PARAMETER      :: N9PFxl         = 265
   INTEGER(IntKi), PARAMETER      :: N9PFyl         = 266
   INTEGER(IntKi), PARAMETER      :: N9PFzl         = 267
   INTEGER(IntKi), PARAMETER      :: N1PMxl         = 268
   INTEGER(IntKi), PARAMETER      :: N1PMyl         = 269
   INTEGER(IntKi), PARAMETER      :: N1PMzl         = 270
   INTEGER(IntKi), PARAMETER      :: N2PMxl         = 271
   INTEGER(IntKi), PARAMETER      :: N2PMyl         = 272
   INTEGER(IntKi), PARAMETER      :: N2PMzl         = 273
   INTEGER(IntKi), PARAMETER      :: N3PMxl         = 274
   INTEGER(IntKi), PARAMETER      :: N3PMyl         = 275
   INTEGER(IntKi), PARAMETER      :: N3PMzl         = 276
   INTEGER(IntKi), PARAMETER      :: N4PMxl         = 277
   INTEGER(IntKi), PARAMETER      :: N4PMyl         = 278
   INTEGER(IntKi), PARAMETER      :: N4PMzl         = 279
   INTEGER(IntKi), PARAMETER      :: N5PMxl         = 280
   INTEGER(IntKi), PARAMETER      :: N5PMyl         = 281
   INTEGER(IntKi), PARAMETER      :: N5PMzl         = 282
   INTEGER(IntKi), PARAMETER      :: N6PMxl         = 283
   INTEGER(IntKi), PARAMETER      :: N6PMyl         = 284
   INTEGER(IntKi), PARAMETER      :: N6PMzl         = 285
   INTEGER(IntKi), PARAMETER      :: N7PMxl         = 286
   INTEGER(IntKi), PARAMETER      :: N7PMyl         = 287
   INTEGER(IntKi), PARAMETER      :: N7PMzl         = 288
   INTEGER(IntKi), PARAMETER      :: N8PMxl         = 289
   INTEGER(IntKi), PARAMETER      :: N8PMyl         = 290
   INTEGER(IntKi), PARAMETER      :: N8PMzl         = 291
   INTEGER(IntKi), PARAMETER      :: N9PMxl         = 292
   INTEGER(IntKi), PARAMETER      :: N9PMyl         = 293
   INTEGER(IntKi), PARAMETER      :: N9PMzl         = 294
   INTEGER(IntKi), PARAMETER      :: N1DFxl         = 295
   INTEGER(IntKi), PARAMETER      :: N1DFyl         = 296
   INTEGER(IntKi), PARAMETER      :: N1DFzl         = 297
   INTEGER(IntKi), PARAMETER      :: N2DFxl         = 298
   INTEGER(IntKi), PARAMETER      :: N2DFyl         = 299
   INTEGER(IntKi), PARAMETER      :: N2DFzl         = 300
   INTEGER(IntKi), PARAMETER      :: N3DFxl         = 301
   INTEGER(IntKi), PARAMETER      :: N3DFyl         = 302
   INTEGER(IntKi), PARAMETER      :: N3DFzl         = 303
   INTEGER(IntKi), PARAMETER      :: N4DFxl         = 304
   INTEGER(IntKi), PARAMETER      :: N4DFyl         = 305
   INTEGER(IntKi), PARAMETER      :: N4DFzl         = 306
   INTEGER(IntKi), PARAMETER      :: N5DFxl         = 307
   INTEGER(IntKi), PARAMETER      :: N5DFyl         = 308
   INTEGER(IntKi), PARAMETER      :: N5DFzl         = 309
   INTEGER(IntKi), PARAMETER      :: N6DFxl         = 310
   INTEGER(IntKi), PARAMETER      :: N6DFyl         = 311
   INTEGER(IntKi), PARAMETER      :: N6DFzl         = 312
   INTEGER(IntKi), PARAMETER      :: N7DFxl         = 313
   INTEGER(IntKi), PARAMETER      :: N7DFyl         = 314
   INTEGER(IntKi), PARAMETER      :: N7DFzl         = 315
   INTEGER(IntKi), PARAMETER      :: N8DFxl         = 316
   INTEGER(IntKi), PARAMETER      :: N8DFyl         = 317
   INTEGER(IntKi), PARAMETER      :: N8DFzl         = 318
   INTEGER(IntKi), PARAMETER      :: N9DFxl         = 319
   INTEGER(IntKi), PARAMETER      :: N9DFyl         = 320
   INTEGER(IntKi), PARAMETER      :: N9DFzl         = 321
   INTEGER(IntKi), PARAMETER      :: N1DMxl         = 322
   INTEGER(IntKi), PARAMETER      :: N1DMyl         = 323
   INTEGER(IntKi), PARAMETER      :: N1DMzl         = 324
   INTEGER(IntKi), PARAMETER      :: N2DMxl         = 325
   INTEGER(IntKi), PARAMETER      :: N2DMyl         = 326
   INTEGER(IntKi), PARAMETER      :: N2DMzl         = 327
   INTEGER(IntKi), PARAMETER      :: N3DMxl         = 328
   INTEGER(IntKi), PARAMETER      :: N3DMyl         = 329
   INTEGER(IntKi), PARAMETER      :: N3DMzl         = 330
   INTEGER(IntKi), PARAMETER      :: N4DMxl         = 331
   INTEGER(IntKi), PARAMETER      :: N4DMyl         = 332
   INTEGER(IntKi), PARAMETER      :: N4DMzl         = 333
   INTEGER(IntKi), PARAMETER      :: N5DMxl         = 334
   INTEGER(IntKi), PARAMETER      :: N5DMyl         = 335
   INTEGER(IntKi), PARAMETER      :: N5DMzl         = 336
   INTEGER(IntKi), PARAMETER      :: N6DMxl         = 337
   INTEGER(IntKi), PARAMETER      :: N6DMyl         = 338
   INTEGER(IntKi), PARAMETER      :: N6DMzl         = 339
   INTEGER(IntKi), PARAMETER      :: N7DMxl         = 340
   INTEGER(IntKi), PARAMETER      :: N7DMyl         = 341
   INTEGER(IntKi), PARAMETER      :: N7DMzl         = 342
   INTEGER(IntKi), PARAMETER      :: N8DMxl         = 343
   INTEGER(IntKi), PARAMETER      :: N8DMyl         = 344
   INTEGER(IntKi), PARAMETER      :: N8DMzl         = 345
   INTEGER(IntKi), PARAMETER      :: N9DMxl         = 346
   INTEGER(IntKi), PARAMETER      :: N9DMyl         = 347
   INTEGER(IntKi), PARAMETER      :: N9DMzl         = 348


     ! Applied loads mapped to root (includes distributed and point):

   INTEGER(IntKi), PARAMETER      :: RootAppliedFxr = 349
   INTEGER(IntKi), PARAMETER      :: RootAppliedFyr = 350
   INTEGER(IntKi), PARAMETER      :: RootAppliedFzr = 351
   INTEGER(IntKi), PARAMETER      :: RootAppliedMxr = 352
   INTEGER(IntKi), PARAMETER      :: RootAppliedMyr = 353
   INTEGER(IntKi), PARAMETER      :: RootAppliedMzr = 354
   INTEGER(IntKi), PARAMETER      :: RootAppliedFxg = 355
   INTEGER(IntKi), PARAMETER      :: RootAppliedFyg = 356
   INTEGER(IntKi), PARAMETER      :: RootAppliedFzg = 357
   INTEGER(IntKi), PARAMETER      :: RootAppliedMxg = 358
   INTEGER(IntKi), PARAMETER      :: RootAppliedMyg = 359
   INTEGER(IntKi), PARAMETER      :: RootAppliedMzg = 360


     ! Pitch Actuator:

   INTEGER(IntKi), PARAMETER      :: PAngInp        = 361
   INTEGER(IntKi), PARAMETER      :: PAngAct        = 362
   INTEGER(IntKi), PARAMETER      :: PRatAct        = 363
   INTEGER(IntKi), PARAMETER      :: PAccAct        = 364


     ! The maximum number of output channels which can be output by the code.
   INTEGER(IntKi), PARAMETER      :: MaxOutPts      = 364

!End of code generated by Matlab script
! ===================================================================================================

   INTEGER,  PARAMETER          :: NFl(9, 3) = RESHAPE( (/ &      ! Sectional force resultants expressed in l
                                     N1Fxl,N2Fxl,N3Fxl,N4Fxl,N5Fxl,N6Fxl,N7Fxl,N8Fxl,N9Fxl, &
                                     N1Fyl,N2Fyl,N3Fyl,N4Fyl,N5Fyl,N6Fyl,N7Fyl,N8Fyl,N9Fyl, &
                                     N1Fzl,N2Fzl,N3Fzl,N4Fzl,N5Fzl,N6Fzl,N7Fzl,N8Fzl,N9Fzl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NMl(9, 3) = RESHAPE( (/ &      ! Sectional moment resultants expressed in l
                                     N1Mxl,N2Mxl,N3Mxl,N4Mxl,N5Mxl,N6Mxl,N7Mxl,N8Mxl,N9Mxl, &
                                     N1Myl,N2Myl,N3Myl,N4Myl,N5Myl,N6Myl,N7Myl,N8Myl,N9Myl, &
                                     N1Mzl,N2Mzl,N3Mzl,N4Mzl,N5Mzl,N6Mzl,N7Mzl,N8Mzl,N9Mzl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NTDr(9, 3) = RESHAPE( (/ &      ! Sectional translational deflections expressed in r
                                     N1TDxr,N2TDxr,N3TDxr,N4TDxr,N5TDxr,N6TDxr,N7TDxr,N8TDxr,N9TDxr, &
                                     N1TDyr,N2TDyr,N3TDyr,N4TDyr,N5TDyr,N6TDyr,N7TDyr,N8TDyr,N9TDyr, &
                                     N1TDzr,N2TDzr,N3TDzr,N4TDzr,N5TDzr,N6TDzr,N7TDzr,N8TDzr,N9TDzr  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NRDr(9, 3) = RESHAPE( (/ &      ! Sectional angular/rotational deflections expressed in r
                                     N1RDxr,N2RDxr,N3RDxr,N4RDxr,N5RDxr,N6RDxr,N7RDxr,N8RDxr,N9RDxr, &
                                     N1RDyr,N2RDyr,N3RDyr,N4RDyr,N5RDyr,N6RDyr,N7RDyr,N8RDyr,N9RDyr, &
                                     N1RDzr,N2RDzr,N3RDzr,N4RDzr,N5RDzr,N6RDzr,N7RDzr,N8RDzr,N9RDzr  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NTVg(9, 3) = RESHAPE( (/ &      ! Sectional translational velocities expressed in g
                                     N1TVXg,N2TVXg,N3TVXg,N4TVXg,N5TVXg,N6TVXg,N7TVXg,N8TVXg,N9TVXg, &
                                     N1TVYg,N2TVYg,N3TVYg,N4TVYg,N5TVYg,N6TVYg,N7TVYg,N8TVYg,N9TVYg, &
                                     N1TVZg,N2TVZg,N3TVZg,N4TVZg,N5TVZg,N6TVZg,N7TVZg,N8TVZg,N9TVZg  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NRVg(9, 3) = RESHAPE( (/ &      ! Sectional angular/rotational velocities expressed in g
                                     N1RVXg,N2RVXg,N3RVXg,N4RVXg,N5RVXg,N6RVXg,N7RVXg,N8RVXg,N9RVXg, &
                                     N1RVYg,N2RVYg,N3RVYg,N4RVYg,N5RVYg,N6RVYg,N7RVYg,N8RVYg,N9RVYg, &
                                     N1RVZg,N2RVZg,N3RVZg,N4RVZg,N5RVZg,N6RVZg,N7RVZg,N8RVZg,N9RVZg  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NTAl(9, 3) = RESHAPE( (/ &      ! Sectional translational accelerations expressed in l
                                     N1TAXl,N2TAXl,N3TAXl,N4TAXl,N5TAXl,N6TAXl,N7TAXl,N8TAXl,N9TAXl, &
                                     N1TAYl,N2TAYl,N3TAYl,N4TAYl,N5TAYl,N6TAYl,N7TAYl,N8TAYl,N9TAYl, &
                                     N1TAZl,N2TAZl,N3TAZl,N4TAZl,N5TAZl,N6TAZl,N7TAZl,N8TAZl,N9TAZl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NRAl(9, 3) = RESHAPE( (/ &      ! Sectional angular/rotational accelerations expressed in l
                                     N1RAXl,N2RAXl,N3RAXl,N4RAXl,N5RAXl,N6RAXl,N7RAXl,N8RAXl,N9RAXl, &
                                     N1RAYl,N2RAYl,N3RAYl,N4RAYl,N5RAYl,N6RAYl,N7RAYl,N8RAYl,N9RAYl, &
                                     N1RAZl,N2RAZl,N3RAZl,N4RAZl,N5RAZl,N6RAZl,N7RAZl,N8RAZl,N9RAZl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NPFl(9, 3) = RESHAPE( (/ &      ! Applied point forces expressed in l
                                     N1PFxl,N2PFxl,N3PFxl,N4PFxl,N5PFxl,N6PFxl,N7PFxl,N8PFxl,N9PFxl, &
                                     N1PFyl,N2PFyl,N3PFyl,N4PFyl,N5PFyl,N6PFyl,N7PFyl,N8PFyl,N9PFyl, &
                                     N1PFzl,N2PFzl,N3PFzl,N4PFzl,N5PFzl,N6PFzl,N7PFzl,N8PFzl,N9PFzl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NPMl(9, 3) = RESHAPE( (/ &      ! Applied point moments expressed in l
                                     N1PMxl,N2PMxl,N3PMxl,N4PMxl,N5PMxl,N6PMxl,N7PMxl,N8PMxl,N9PMxl, &
                                     N1PMyl,N2PMyl,N3PMyl,N4PMyl,N5PMyl,N6PMyl,N7PMyl,N8PMyl,N9PMyl, &
                                     N1PMzl,N2PMzl,N3PMzl,N4PMzl,N5PMzl,N6PMzl,N7PMzl,N8PMzl,N9PMzl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NDFl(9, 3) = RESHAPE( (/ &      ! Applied distributed forces expressed in l
                                     N1DFxl,N2DFxl,N3DFxl,N4DFxl,N5DFxl,N6DFxl,N7DFxl,N8DFxl,N9DFxl, &
                                     N1DFyl,N2DFyl,N3DFyl,N4DFyl,N5DFyl,N6DFyl,N7DFyl,N8DFyl,N9DFyl, &
                                     N1DFzl,N2DFzl,N3DFzl,N4DFzl,N5DFzl,N6DFzl,N7DFzl,N8DFzl,N9DFzl  &
                                   /), (/9, 3/) )
   INTEGER,  PARAMETER          :: NDMl(9, 3) = RESHAPE( (/ &      ! Applied distributed moments expressed in l
                                     N1DMxl,N2DMxl,N3DMxl,N4DMxl,N5DMxl,N6DMxl,N7DMxl,N8DMxl,N9DMxl, &
                                     N1DMyl,N2DMyl,N3DMyl,N4DMyl,N5DMyl,N6DMyl,N7DMyl,N8DMyl,N9DMyl, &
                                     N1DMzl,N2DMzl,N3DMzl,N4DMzl,N5DMzl,N6DMzl,N7DMzl,N8DMzl,N9DMzl  &
                                   /), (/9, 3/) )


CONTAINS
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE BD_ReadInput(InputFileName,InputFileData,OutFileRoot, Default_DT,ErrStat,ErrMsg)

   ! Passed Variables:
   CHARACTER(*),                 INTENT(IN   )  :: InputFileName    !< Name of the input file
   CHARACTER(*),                 INTENT(IN   )  :: OutFileRoot      !< Name of the input file
   REAL(DbKi),                   INTENT(IN   )  :: Default_DT       !< The default DT (from glue code)
   TYPE(BD_InputFile),           INTENT(  OUT)  :: InputFileData    !< Data stored in the module's input file
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat          !< The error status code
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg           !< The error message, if an error occurred


   ! Local variables:
   INTEGER(IntKi)                               :: UnEcho
   INTEGER(IntKi)                               :: ErrStat2
   CHARACTER(ErrMsgLen)                         :: ErrMsg2
   character(*), parameter                      :: RoutineName = 'BD_ReadInput'


   ErrStat = ErrID_None
   ErrMsg = ''
   UnEcho = -1

   InputFileData%DTBeam = Default_DT
   CALL BD_ReadPrimaryFile(InputFileName,InputFileData,OutFileRoot,UnEcho,ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
      IF(ErrStat >= AbortErrLev) THEN
         IF (UnEcho > 0) CLOSE (UnEcho)
         RETURN
      END IF

   CALL BD_ReadBladeFile(InputFileData%BldFile,InputFileData%InpBl,UnEcho,ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
      IF(ErrStat >= AbortErrLev) THEN
         IF (UnEcho > 0) CLOSE (UnEcho)
         RETURN
      END IF

END SUBROUTINE BD_ReadInput
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine reads in the primary BeamDyn input file and places the values it reads
!! in the InputFileData structure.
!!   It opens an echo file if requested and returns the (still-open) echo file to the
!!     calling routine.
SUBROUTINE BD_ReadPrimaryFile(InputFile,InputFileData,OutFileRoot,UnEc,ErrStat,ErrMsg)

   ! Passed variables
   INTEGER(IntKi),               INTENT(  OUT) :: UnEc
   INTEGER(IntKi),               INTENT(  OUT) :: ErrStat
   CHARACTER(*),                 INTENT(IN   ) :: InputFile
   CHARACTER(*),                 INTENT(IN   ) :: OutFileRoot
   CHARACTER(*),                 INTENT(  OUT) :: ErrMsg

   TYPE(BD_InputFile),           INTENT(INOUT) :: InputFileData

   ! Local variables:
   INTEGER(IntKi)               :: UnIn                         ! Unit number for reading file
   INTEGER(IntKi)               :: ErrStat2                     ! Temporary Error status
   LOGICAL                      :: Echo                         ! Determines if an echo file should be written
   INTEGER(IntKi)               :: IOS                          ! Temporary Error status
   CHARACTER(ErrMsgLen)         :: ErrMsg2                      ! Temporary Error message
   CHARACTER(ErrMsgLen)         :: ErrMsg_NoBldNdOuts           ! Temporary Error message
   character(*), parameter      :: RoutineName = 'BD_ReadPrimaryFile'

   CHARACTER(1024)              :: PriPath                      ! Path name of the primary file
   CHARACTER(1024)              :: FTitle                       ! "File Title": the 2nd line of the input file, which contains a description of its contents
   CHARACTER(200)               :: Line                         ! Temporary storage of a line from the input

   INTEGER(IntKi)               :: i
   INTEGER(IntKi)               :: j
   INTEGER(IntKi)               :: temp_int
   REAL(BDKi)                   :: tmpReAry(4)

   ! Initialize some variables:
   ErrStat = ErrID_None
   ErrMsg  = ""
   Echo = .FALSE.
   UnEc = -1
   CALL GetPath( InputFile, PriPath )     ! Input files will be relative to the path where the primary input file is located.

   CALL AllocAry( InputFileData%OutList, MaxOutPts, "Outlist", ErrStat2, ErrMsg2 )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )


      ! Allocate array for holding the list of node outputs
   CALL AllocAry( InputFileData%BldNd_OutList, BldNd_MaxOutPts, "BldNd_Outlist", ErrStat2, ErrMsg2 )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )


   !$OMP critical(fileopen_critical)
   CALL GetNewUnit(UnIn,ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL OpenFInpFile(UnIn,InputFile,ErrStat2,ErrMsg2)
   !$OMP end critical(fileopen_critical)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      IF (ErrStat >= AbortErrLev) then
         call cleanup()
         RETURN
      end if



   I = 1 !set the number of times we've read the file
   DO
   !-------------------------- HEADER ---------------------------------------------
   CALL ReadCom(UnIn,InputFile,'File Header: Module Version (line 1)',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadStr(UnIn,InputFile,FTitle,'FTitle','File Header: File Description (line 2)',ErrStat2, ErrMsg2, UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if


   !---------------------- SIMULATION CONTROL --------------------------------------
   CALL ReadCom(UnIn,InputFile,'Section Header: Simulation Control',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,InputFile,Echo,'Echo','Echo switch',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if


      IF (.NOT. Echo .OR. I > 1) EXIT !exit this loop

         ! Otherwise, open the echo file, then rewind the input file and echo everything we've read

      I = I + 1         ! make sure we do this only once (increment counter that says how many times we've read this file)

      CALL OpenEcho ( UnEc, TRIM(OutFileRoot)//'.ech', ErrStat2, ErrMsg2, BeamDyn_Ver )
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         if (ErrStat >= AbortErrLev) then
            call cleanup()
            return
         end if

      IF ( UnEc > 0 )  WRITE (UnEc,'(/,A,/)')  'Data from '//TRIM(BeamDyn_Ver%Name)//' primary input file "'//TRIM( InputFile )//'":'

      REWIND( UnIn, IOSTAT=ErrStat2 )
         IF (ErrStat2 /= 0_IntKi ) THEN
            CALL SetErrStat( ErrID_Fatal, 'Error rewinding file "'//TRIM(InputFile)//'".', ErrStat, ErrMsg, RoutineName )
            call cleanup()
            RETURN
         END IF

   END DO

   CALL ReadVar(UnIn,InputFile,InputFileData%QuasiStaticInit,"QuasiStaticInit", "QuasiStaticInit",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,InputFile,InputFileData%rhoinf,"rhoinf", "Coefficient for GA2",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,InputFile,InputFileData%quadrature,"quadrature", "Quadrature type",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "refine", "Refinement parameter for trapezoidal quadrature {or default} ", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
          InputFileData%refine = 1
      ELSE ! If it's not "default", read this variable
         READ( Line, *, IOSTAT=IOS) InputFileData%refine
            CALL CheckIOS ( IOS, InputFile, 'refine', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "n_fact", "Factorization frequency {or default}", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
          InputFileData%n_fact = 5
      ELSE ! If it's not "default", read this variable
         READ( Line, *, IOSTAT=IOS) InputFileData%n_fact
            CALL CheckIOS ( IOS, InputFile, 'n_fact', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "DTBeam", "Time interval for BeamDyn  calculations {or default} (s)", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) /= 1 ) THEN ! If it's not "default", read this variable; otherwise use the value already stored in InputFileData%DTBeam
         READ( Line, *, IOSTAT=IOS) InputFileData%DTBeam
            CALL CheckIOS ( IOS, InputFile, 'DTBeam', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "load_retries", "Tolerance for stopping criterion", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
          InputFileData%load_retries = 20
      ELSE ! If it's not "default", read this variable;
         READ( Line, *, IOSTAT=IOS) InputFileData%load_retries
            CALL CheckIOS ( IOS, InputFile, 'load_retries', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "NRMax", "Max number of interations in Newton-Raphson algorithm", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
          InputFileData%NRMax = 10
      ELSE ! If it's not "default", read this variable;
         READ( Line, *, IOSTAT=IOS) InputFileData%NRMax
            CALL CheckIOS ( IOS, InputFile, 'NRMax', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar( UnIn, InputFile, Line, "stop_tol", "Tolerance for stopping criterion", ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if
      CALL Conv2UC( Line )
      IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
          InputFileData%stop_tol = 1.0D-05
      ELSE ! If it's not "default", read this variable;
         READ( Line, *, IOSTAT=IOS) InputFileData%stop_tol
            CALL CheckIOS ( IOS, InputFile, 'stop_tol', NumType, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      END IF

   Line = ""
   CALL ReadVar(UnIn, InputFile, Line, 'tngt_stf_fd','finite difference for tangent stiffness flag', ErrStat2, ErrMsg2, UnEc)
   CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   if (ErrStat >= AbortErrLev) then
      call cleanup()
      return
   end if
   CALL Conv2UC( Line )
   IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
       InputFileData%tngt_stf_fd = .FALSE.
   ELSE ! If it's not "default", read this variable;
       READ( Line, *, IOSTAT=IOS) InputFileData%tngt_stf_fd
       CALL CheckIOS ( IOS, InputFile, 'tngt_stf_fd', NumType, ErrStat2, ErrMsg2 )
       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   END IF
   if (InputFileData%tngt_stf_fd) CALL WrScr( 'Using finite difference to compute tangent stiffness matrix'//NewLine )
   !   ! RelStates - Define states relative to root motion during linearization? (flag) [used only when linearizing]
   !CALL ReadVar(UnIn,InputFile,InputFileData%RelStates,"RelStates", "Define states relative to root motion during linearization? (flag) [used only when linearizing]",ErrStat2,ErrMsg2,UnEc)
   !   CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   InputFileData%RelStates = .false.  ! this doesn't seem to be needed anymore (and I think there is a problem with using it in MBC3)

   Line = ""
   CALL ReadVar(UnIn, InputFile, Line, 'tngt_stf_comp','compare tangent stiffness using finite difference flag', ErrStat2, ErrMsg2, UnEc)
   CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   if (ErrStat >= AbortErrLev) then
      call cleanup()
      return
   end if
   CALL Conv2UC( Line )
   IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
       InputFileData%tngt_stf_comp = .FALSE.
   ELSE ! If it's not "default", read this variable;
       READ( Line, *, IOSTAT=IOS) InputFileData%tngt_stf_comp
       CALL CheckIOS ( IOS, InputFile, 'tngt_stf_comp', NumType, ErrStat2, ErrMsg2 )
       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   END IF
   if (InputFileData%tngt_stf_comp) CALL WrScr( 'WARNING: tngt_stf_comp set to true. Output will be verbose'//NewLine )

   Line = ""
   CALL ReadVar(UnIn, InputFile, Line, 'tngt_stf_pert','perturbation size for finite differenced tangent stiffness', ErrStat2, ErrMsg2, UnEc)
   CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   if (ErrStat >= AbortErrLev) then
      call cleanup()
      return
   end if
   CALL Conv2UC( Line )
   IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
       InputFileData%tngt_stf_pert = 1.0D-06
   ELSE ! If it's not "default", read this variable;
       READ( Line, *, IOSTAT=IOS) InputFileData%tngt_stf_pert
       CALL CheckIOS ( IOS, InputFile, 'tngt_stf_pert', NumType, ErrStat2, ErrMsg2 )
       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   END IF

   Line = ""
   CALL ReadVar(UnIn, InputFile, Line, 'tngt_stf_difftol','tolerance for difference between analytical and fd tangent stiffness', ErrStat2, ErrMsg2, UnEc)
   CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   if (ErrStat >= AbortErrLev) then
      call cleanup()
      return
   end if
   CALL Conv2UC( Line )
   IF ( INDEX(Line, "DEFAULT" ) .EQ. 1) THEN
       InputFileData%tngt_stf_difftol = 1.0D-01
   ELSE ! If it's not "default", read this variable;
       READ( Line, *, IOSTAT=IOS) InputFileData%tngt_stf_difftol
       CALL CheckIOS ( IOS, InputFile, 'tngt_stf_difftol', NumType, ErrStat2, ErrMsg2 )
       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   END IF

   ! RotStates - Orient states in the rotating frame during linearization? (flag) [used only when linearizing]
   CALL ReadVar(UnIn,InputFile,InputFileData%RotStates,"RotStates", "Orient states in the rotating frame during linearization? (flag) [used only when linearizing]",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if

   ! return on error at end of section
   if (ErrStat >= AbortErrLev) then
       call cleanup()
       return
   end if

   !---------------------- GEOMETRY PARAMETER --------------------------------------
   CALL ReadCom(UnIn,InputFile,'Section Header: Geometry Parameter',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,InputFile,InputFileData%member_total,"member_total", "Total number of member",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,InputFile,InputFileData%kp_total,"kp_total", "Total number of key point",ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

         ! check these values before allocating the arrays below
      if (InputFileData%member_total<1 .or. InputFileData%kp_total<1) then
         CALL SetErrStat( ErrID_Fatal, "member_total and kp_total must be positive numbers", ErrStat, ErrMsg, RoutineName )
         call cleanup()
         return
      end if

      CALL AllocAry(InputFileData%kp_member,InputFileData%member_total,'Number of key point in each member',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      CALL AllocAry(InputFileData%kp_coordinate,InputFileData%kp_total,4,'Key point coordinates input array',ErrStat2,ErrMsg2)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         if (ErrStat >= AbortErrLev) then
            call cleanup()
            return
         end if

   InputFileData%kp_member = 0
   temp_int = 0
   DO i=1,InputFileData%member_total
      ! bjj: we cannot read j, InputFileData%kp_member(j) because j could be outside the valid range:
      READ(UnIn,*,IOSTAT=IOS) j,temp_int
         CALL CheckIOS ( IOS, InputFile, 'Member number; Number of key points in this member', NumType, ErrStat2, ErrMsg2 )
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         if (ErrStat >= AbortErrLev) then
            call cleanup()
            return
         end if

         if (j>0 .and. j<= InputFileData%member_total) then
            InputFileData%kp_member(j) = temp_int
         end if

         !bjj: what happens if the user enters a member number multiple times? or if not each member is accounted for?
         !    I'm going to force them to be entered 1-InputFileData%kp_member to avoid this issue for now.
         if (j /= i) then
            call SetErrStat(ErrID_Warn, "Member numbers must be entered in monotonicly increasing order, starting with 1.",ErrStat,ErrMsg,RoutineName)
         end if

   ENDDO


   CALL ReadCom(UnIn,InputFile,'key point x,y,z locations and initial twist angles',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadCom(UnIn,InputFile,'key point and initial twist units',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   !InputFileData%kp_coordinate = 0.0_BDKi
   DO i=1,InputFileData%kp_total
       CALL ReadAry( UnIn, InputFile, TmpReAry, 4, 'kp_coordinate', 'Key point coordinates and initial twist', ErrStat2, ErrMsg2, UnEc )
          CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
       InputFileData%kp_coordinate(i,1) =  TmpReAry(1) ! x
       InputFileData%kp_coordinate(i,2) =  TmpReAry(2) ! y
       InputFileData%kp_coordinate(i,3) =  TmpReAry(3) ! z
       InputFileData%kp_coordinate(i,4) = -TmpReAry(4) ! initial twist
   ENDDO


   !---------------------- MESH PARAMETER -----------------------------------------
   CALL ReadCom(UnIn,InputFile,'Section Header: Mesh Parameter',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL ReadVar(UnIn,InputFile,InputFileData%order_elem,"order_elem","Order of basis function",&
                ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   !---------------------- BEAM SECTIONAL PARAMETER ----------------------------------------
   CALL ReadCom(UnIn,InputFile,'Section Header: Beam Material Parameter',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL ReadVar ( UnIn, InputFile, InputFileData%BldFile, 'BldFile', 'Name of the file containing properties for beam', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      IF ( PathIsRelative( InputFileData%BldFile ) ) InputFileData%BldFile = TRIM(PriPath)//TRIM(InputFileData%BldFile)


   !---------------------- PITCH ACTUATOR PARAMETER ----------------------------------------
   CALL ReadCom(UnIn,InputFile,'Section Header: Pitch Actuator Parameter',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar ( UnIn, InputFile, InputFileData%UsePitchAct, 'UsePitchAct', 'Whether a pitch actuator should be used (flag)', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar ( UnIn, InputFile, InputFileData%PitchJ, 'PitchJ', 'Pitch actuator inertia (kg-m^2)', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar ( UnIn, InputFile, InputFileData%PitchK, 'PitchK', 'Pitch actuator stiffness (kg-m^2/s^2)', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar ( UnIn, InputFile, InputFileData%PitchC, 'PitchC', 'Pitch actuator damping (kg-m^2/s)', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )


   !----------- OUTPUTS  -----------------------------------------------------------
   CALL ReadCom( UnIn, InputFile, 'Section Header: Outputs', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      ! SumPrint - Print summary data to <RootName>.BD.sum (flag):
   CALL ReadVar( UnIn, InputFile, InputFileData%SumPrint, "SumPrint", "Print summary data to <RootName>.BD.sum (flag)", ErrStat2, ErrMsg2, UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   call ReadVar( UnIn, InputFile, InputFileData%OutFmt, "OutFmt", "Format used for text tabular output (except time).  Resulting field should be 10 characters. (-)",ErrStat2, ErrMsg2, UnEc)
      call setErrStat( errStat2, ErrMsg2 , errStat, ErrMsg , RoutineName )

      ! NNodeOuts - Number of node outputs [0 - 9] (-):
   CALL ReadVar( UnIn, InputFile, InputFileData%NNodeOuts, "NNodeOuts", "Number of node outputs [0 - 9] (-)", ErrStat2, ErrMsg2, UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      IF ( InputFileData%NNodeOuts > SIZE(InputFileData%OutNd) ) THEN
         CALL SetErrStat( ErrID_Warn, ' Warning: number of output nodes exceeds '//&
                           TRIM(Num2LStr(SIZE(InputFileData%OutNd))) //'.', ErrStat, ErrMsg, RoutineName )
         InputFileData%NNodeOuts = SIZE(InputFileData%OutNd)
      END IF

      ! OutNd - Nodes whose values will be output (-):
   CALL ReadAry( UnIn, InputFile, InputFileData%OutNd, InputFileData%NNodeOuts, "OutNd", "Nodes whose values will be output (-)", ErrStat2, ErrMsg2, UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )


      ! Return on error at end of section
   IF ( ErrStat >= AbortErrLev ) THEN
      CALL Cleanup()
      RETURN
   END IF

   !----------- OUTLIST  -----------------------------------------------------------
   CALL ReadCom( UnIn, InputFile, 'Section Header: OutList', ErrStat2, ErrMsg2, UnEc )
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      ! OutList - List of user-requested output channels (-):
   CALL ReadOutputList ( UnIn, InputFile, InputFileData%OutList, InputFileData%NumOuts, 'OutList', "List of user-requested output channels", ErrStat2, ErrMsg2, UnEc  )     ! Routine in NWTC Subroutine Library
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   !---------------------- END OF FILE -----------------------------------------



   !----------- OUTLIST  -----------------------------------------------------------
      ! In case there is something ill-formed in the additional nodal outputs section, we will simply ignore it.
   ErrMsg_NoBldNdOuts='Nodal outputs section of BeamDyn input file not found or improperly formatted.'
   InputFileData%BldNd_NumOuts = 0     ! Just in case we don't get an error but have no nodal outputs.


   !----------- OUTLIST for BldNd -----------------------------------------------------------
   CALL ReadCom( UnIn, InputFile, 'Section Header: OutList for Blade node channels', ErrStat2, ErrMsg2, UnEc )
   IF ( ErrStat2 >= AbortErrLev ) THEN
      InputFileData%BldNd_NumOuts = 0
      call wrscr( trim(ErrMsg_NoBldNdOuts) )
      CALL Cleanup()
      RETURN
   ENDIF

      ! Number of blade nodes to output:  will modify this at some point for arrays
      ! TODO: Parse this string into an array of nodes to output at (one idea is to set an array of boolean to T/F for which nodes to output).  At present, we ignore it entirely.
   CALL ReadVar(  UnIn, InputFile, InputFileData%BldNd_BlOutNd_Str, 'BldNd_BlOutNd_Str', 'Which nodes to output node data on.'//TRIM(Num2Lstr(I)), ErrStat2, ErrMsg2, UnEc )
   IF ( ErrStat2 >= AbortErrLev ) THEN
      InputFileData%BldNd_NumOuts = 0
      call wrscr( trim(ErrMsg_NoBldNdOuts) )
      CALL Cleanup()
      RETURN
   ENDIF


      ! Section header for outlist
   CALL ReadCom( UnIn, InputFile, 'Section Header: OutList', ErrStat2, ErrMsg2, UnEc )
   IF ( ErrStat2 >= AbortErrLev ) THEN
      InputFileData%BldNd_NumOuts = 0
      CALL SetErrStat( ErrID_Warn, ErrMsg_NoBldNdOuts, ErrStat, ErrMsg, RoutineName )
      CALL Cleanup()
      RETURN
   ENDIF


      ! OutList - List of user-requested output channels at each node(-):
   CALL ReadOutputList ( UnIn, InputFile, InputFileData%BldNd_OutList, InputFileData%BldNd_NumOuts, 'BldNd_OutList', "List of user-requested output channels", ErrStat2, ErrMsg2, UnEc  )     ! Routine in NWTC Subroutine Library
   IF ( ErrStat2 >= AbortErrLev .and. InputFileData%BldNd_NumOuts < 1) THEN
      InputFileData%BldNd_NumOuts = 0
      call wrscr( trim(ErrMsg_NoBldNdOuts) )
      CALL Cleanup()
      RETURN
   ELSE
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   ENDIF

   !---------------------- END OF FILE -----------------------------------------


   call cleanup()
   return

contains
   subroutine cleanup()
      if (UnIn > 0) close(UnIn)
      ! don't close the echo file here
   end subroutine cleanup
END SUBROUTINE BD_ReadPrimaryFile
!----------------------------------------------------------------------------------------------------------------------------------
SUBROUTINE BD_ReadBladeFile(BldFile,BladeInputFileData,UnEc,ErrStat,ErrMsg)

   ! Passed variables:
   TYPE(BladeInputData), INTENT(  OUT):: BladeInputFileData
   CHARACTER(*),         INTENT(IN   ):: BldFile
   INTEGER(IntKi),       INTENT(IN   ):: UnEc
   INTEGER(IntKi),       INTENT(  OUT):: ErrStat                             !< Error status
   CHARACTER(*),         INTENT(  OUT):: ErrMsg                              !< Error message

   ! Local variables:
   INTEGER(IntKi)             :: UnIn                                            ! Unit number for reading file
   INTEGER(IntKi)             :: ErrStat2                                        ! Temporary Error status
   CHARACTER(ErrMsgLen)       :: ErrMsg2                                         ! Temporary Err msg
   character(*), parameter    :: RoutineName = 'BD_ReadBladeFile'
   INTEGER(IntKi)             :: i
   INTEGER(IntKi)             :: j

   REAL(BDKi)                 :: temp66(6,6)

   ErrStat = ErrID_None
   ErrMsg  = ""

   !$OMP critical(fileopen_critical)
   CALL GetNewUnit(UnIn,ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL OpenFInpFile (UnIn,BldFile,ErrStat2,ErrMsg2)
   !$OMP end critical(fileopen_critical)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
         if (ErrStat >= AbortErrLev) then
            return
         end if

   !  -------------- HEADER -------------------------------------------------------
   ! Skip the header.
   CALL ReadCom(UnIn,BldFile,'unused beam file header line 1',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL ReadCom(UnIn,BldFile,'unused beam file header line 2',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   !  -------------- BLADE PARAMETER-----------------------------------------------
   CALL ReadCom(UnIn,BldFile,'beam parameters',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadVar(UnIn,BldFile,BladeInputFileData%station_total,'station_total','Number of blade input stations',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL AllocAry(BladeInputFileData%stiff0,6,6,BladeInputFileData%station_total,'Cross-sectional 6 by 6 stiffness matrix',ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL AllocAry(BladeInputFileData%mass0,6,6,BladeInputFileData%station_total,'Cross-sectional 6 by 6 mass matrix',ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL AllocAry(BladeInputFileData%station_eta,BladeInputFileData%station_total,'Station eta array',ErrStat2,ErrMsg2)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      !after allocating these arrays, we'll make sure it was successful:
   if (ErrStat >= AbortErrLev) then
      call cleanup()
      return
   end if

!FIXME: any reason not to use a logical for this?
   CALL ReadVar(UnIn,BldFile,BladeInputFileData%damp_flag,'damp_flag','Damping flag',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   !  -------------- DAMPING PARAMETER-----------------------------------------------
   CALL ReadCom(UnIn,BldFile,'damping parameters',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL ReadCom(UnIn,BldFile,'mu1 to mu6',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
   CALL ReadCom(UnIn,BldFile,'units',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   CALL ReadAry(UnIn,BldFile,BladeInputFileData%beta,6,'damping coefficient','damping coefficient',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if


!  -------------- DISTRIBUTED PROPERTIES--------------------------------------------
   CALL ReadCom(UnIn,BldFile,'Distributed properties',ErrStat2,ErrMsg2,UnEc)
      CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

   DO i=1,BladeInputFileData%station_total
      if (i > 1) then
         CALL ReadCom(UnIn,BldFile,'blank line after mass matrix',ErrStat2,ErrMsg2,UnEc)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      end if

      CALL ReadVar(UnIn,BldFile,BladeInputFileData%station_eta(i),'station_eta','Station '//trim(num2lstr(i))//' Eta',ErrStat2,ErrMsg2,UnEc)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      DO j=1,6
         CALL ReadAry(UnIn,BldFile,temp66(j,:),6,'siffness_matrix','Blade C/S stiffness matrix',ErrStat2,ErrMsg2,UnEc)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      ENDDO
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if
      BladeInputFileData%stiff0(:,:,i) = temp66
      CALL ReadCom(UnIn,BldFile,'blank line after stiffness matrix',ErrStat2,ErrMsg2,UnEc)
         CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      DO j=1,6
         CALL ReadAry(UnIn,BldFile,temp66(j,:),6,'mass_matrix','Blade C/S mass matrix',ErrStat2,ErrMsg2,UnEc)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
      ENDDO
      if (ErrStat >= AbortErrLev) then
         call cleanup()
         return
      end if
      BladeInputFileData%mass0(:,:,i) = temp66

   ENDDO

   call cleanup()
   return

contains
   !.....................
   subroutine cleanup()
      close(UnIn)
      return
   end subroutine cleanup
END SUBROUTINE BD_ReadBladeFile
!----------------------------------------------------------------------------------------------------------------------------------
!**********************************************************************************************************************************
! NOTE: The following lines of code were generated by a Matlab script called "Write_ChckOutLst.m"
!      using the parameters listed in the "OutListParameters.xlsx" Excel file. Any changes to these 
!      lines should be modified in the Matlab script and/or Excel worksheet as necessary. 
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine checks to see if any requested output channel names (stored in the OutList(:)) are invalid. It returns a 
!! warning if any of the channels are not available outputs from the module.
!!  It assigns the settings for OutParam(:) (i.e, the index, name, and units of the output channels, WriteOutput(:)).
!!  the sign is set to 0 if the channel is invalid.
!! It sets assumes the value p%NumOuts has been set before this routine has been called, and it sets the values of p%OutParam here.
!! 
!! This routine was generated by Write_ChckOutLst.m using the parameters listed in OutListParameters.xlsx at 08-May-2025 15:51:30.
SUBROUTINE SetOutParam(OutList, p, ErrStat, ErrMsg )

   IMPLICIT                        NONE

      ! Passed variables

   CHARACTER(ChanLen),        INTENT(IN)     :: OutList(:)                        !< The list out user-requested outputs
   TYPE(BD_ParameterType),    INTENT(INOUT)  :: p                                 !< The module parameters
   INTEGER(IntKi),            INTENT(OUT)    :: ErrStat                           !< The error status code
   CHARACTER(*),              INTENT(OUT)    :: ErrMsg                            !< The error message, if an error occurred

      ! Local variables

   INTEGER                      :: ErrStat2                                        ! temporary (local) error status
   INTEGER                      :: I                                               ! Generic loop-counting index
!   INTEGER                      :: J                                               ! Generic loop-counting index
   INTEGER                      :: INDX                                            ! Index for valid arrays

   LOGICAL                      :: CheckOutListAgain                               ! Flag used to determine if output parameter starting with "M" is valid (or the negative of another parameter)
   LOGICAL                      :: InvalidOutput(0:MaxOutPts)                      ! This array determines if the output channel is valid for this configuration
   CHARACTER(ChanLen)           :: OutListTmp                                      ! A string to temporarily hold OutList(I)
   CHARACTER(*), PARAMETER      :: RoutineName = "SetOutParam"

   CHARACTER(OutStrLenM1), PARAMETER  :: ValidParamAry(364) =  (/  &   ! This lists the names of the allowed parameters, which must be sorted alphabetically
                               "N1DFXL        ","N1DFYL        ","N1DFZL        ","N1DMXL        ","N1DMYL        ", &
                               "N1DMZL        ","N1FXL         ","N1FYL         ","N1FZL         ","N1MXL         ", &
                               "N1MYL         ","N1MZL         ","N1PFXL        ","N1PFYL        ","N1PFZL        ", &
                               "N1PMXL        ","N1PMYL        ","N1PMZL        ","N1RAXL        ","N1RAYL        ", &
                               "N1RAZL        ","N1RDXR        ","N1RDYR        ","N1RDZR        ","N1RVXG        ", &
                               "N1RVYG        ","N1RVZG        ","N1TAXL        ","N1TAYL        ","N1TAZL        ", &
                               "N1TDXR        ","N1TDYR        ","N1TDZR        ","N1TVXG        ","N1TVYG        ", &
                               "N1TVZG        ","N2DFXL        ","N2DFYL        ","N2DFZL        ","N2DMXL        ", &
                               "N2DMYL        ","N2DMZL        ","N2FXL         ","N2FYL         ","N2FZL         ", &
                               "N2MXL         ","N2MYL         ","N2MZL         ","N2PFXL        ","N2PFYL        ", &
                               "N2PFZL        ","N2PMXL        ","N2PMYL        ","N2PMZL        ","N2RAXL        ", &
                               "N2RAYL        ","N2RAZL        ","N2RDXR        ","N2RDYR        ","N2RDZR        ", &
                               "N2RVXG        ","N2RVYG        ","N2RVZG        ","N2TAXL        ","N2TAYL        ", &
                               "N2TAZL        ","N2TDXR        ","N2TDYR        ","N2TDZR        ","N2TVXG        ", &
                               "N2TVYG        ","N2TVZG        ","N3DFXL        ","N3DFYL        ","N3DFZL        ", &
                               "N3DMXL        ","N3DMYL        ","N3DMZL        ","N3FXL         ","N3FYL         ", &
                               "N3FZL         ","N3MXL         ","N3MYL         ","N3MZL         ","N3PFXL        ", &
                               "N3PFYL        ","N3PFZL        ","N3PMXL        ","N3PMYL        ","N3PMZL        ", &
                               "N3RAXL        ","N3RAYL        ","N3RAZL        ","N3RDXR        ","N3RDYR        ", &
                               "N3RDZR        ","N3RVXG        ","N3RVYG        ","N3RVZG        ","N3TAXL        ", &
                               "N3TAYL        ","N3TAZL        ","N3TDXR        ","N3TDYR        ","N3TDZR        ", &
                               "N3TVXG        ","N3TVYG        ","N3TVZG        ","N4DFXL        ","N4DFYL        ", &
                               "N4DFZL        ","N4DMXL        ","N4DMYL        ","N4DMZL        ","N4FXL         ", &
                               "N4FYL         ","N4FZL         ","N4MXL         ","N4MYL         ","N4MZL         ", &
                               "N4PFXL        ","N4PFYL        ","N4PFZL        ","N4PMXL        ","N4PMYL        ", &
                               "N4PMZL        ","N4RAXL        ","N4RAYL        ","N4RAZL        ","N4RDXR        ", &
                               "N4RDYR        ","N4RDZR        ","N4RVXG        ","N4RVYG        ","N4RVZG        ", &
                               "N4TAXL        ","N4TAYL        ","N4TAZL        ","N4TDXR        ","N4TDYR        ", &
                               "N4TDZR        ","N4TVXG        ","N4TVYG        ","N4TVZG        ","N5DFXL        ", &
                               "N5DFYL        ","N5DFZL        ","N5DMXL        ","N5DMYL        ","N5DMZL        ", &
                               "N5FXL         ","N5FYL         ","N5FZL         ","N5MXL         ","N5MYL         ", &
                               "N5MZL         ","N5PFXL        ","N5PFYL        ","N5PFZL        ","N5PMXL        ", &
                               "N5PMYL        ","N5PMZL        ","N5RAXL        ","N5RAYL        ","N5RAZL        ", &
                               "N5RDXR        ","N5RDYR        ","N5RDZR        ","N5RVXG        ","N5RVYG        ", &
                               "N5RVZG        ","N5TAXL        ","N5TAYL        ","N5TAZL        ","N5TDXR        ", &
                               "N5TDYR        ","N5TDZR        ","N5TVXG        ","N5TVYG        ","N5TVZG        ", &
                               "N6DFXL        ","N6DFYL        ","N6DFZL        ","N6DMXL        ","N6DMYL        ", &
                               "N6DMZL        ","N6FXL         ","N6FYL         ","N6FZL         ","N6MXL         ", &
                               "N6MYL         ","N6MZL         ","N6PFXL        ","N6PFYL        ","N6PFZL        ", &
                               "N6PMXL        ","N6PMYL        ","N6PMZL        ","N6RAXL        ","N6RAYL        ", &
                               "N6RAZL        ","N6RDXR        ","N6RDYR        ","N6RDZR        ","N6RVXG        ", &
                               "N6RVYG        ","N6RVZG        ","N6TAXL        ","N6TAYL        ","N6TAZL        ", &
                               "N6TDXR        ","N6TDYR        ","N6TDZR        ","N6TVXG        ","N6TVYG        ", &
                               "N6TVZG        ","N7DFXL        ","N7DFYL        ","N7DFZL        ","N7DMXL        ", &
                               "N7DMYL        ","N7DMZL        ","N7FXL         ","N7FYL         ","N7FZL         ", &
                               "N7MXL         ","N7MYL         ","N7MZL         ","N7PFXL        ","N7PFYL        ", &
                               "N7PFZL        ","N7PMXL        ","N7PMYL        ","N7PMZL        ","N7RAXL        ", &
                               "N7RAYL        ","N7RAZL        ","N7RDXR        ","N7RDYR        ","N7RDZR        ", &
                               "N7RVXG        ","N7RVYG        ","N7RVZG        ","N7TAXL        ","N7TAYL        ", &
                               "N7TAZL        ","N7TDXR        ","N7TDYR        ","N7TDZR        ","N7TVXG        ", &
                               "N7TVYG        ","N7TVZG        ","N8DFXL        ","N8DFYL        ","N8DFZL        ", &
                               "N8DMXL        ","N8DMYL        ","N8DMZL        ","N8FXL         ","N8FYL         ", &
                               "N8FZL         ","N8MXL         ","N8MYL         ","N8MZL         ","N8PFXL        ", &
                               "N8PFYL        ","N8PFZL        ","N8PMXL        ","N8PMYL        ","N8PMZL        ", &
                               "N8RAXL        ","N8RAYL        ","N8RAZL        ","N8RDXR        ","N8RDYR        ", &
                               "N8RDZR        ","N8RVXG        ","N8RVYG        ","N8RVZG        ","N8TAXL        ", &
                               "N8TAYL        ","N8TAZL        ","N8TDXR        ","N8TDYR        ","N8TDZR        ", &
                               "N8TVXG        ","N8TVYG        ","N8TVZG        ","N9DFXL        ","N9DFYL        ", &
                               "N9DFZL        ","N9DMXL        ","N9DMYL        ","N9DMZL        ","N9FXL         ", &
                               "N9FYL         ","N9FZL         ","N9MXL         ","N9MYL         ","N9MZL         ", &
                               "N9PFXL        ","N9PFYL        ","N9PFZL        ","N9PMXL        ","N9PMYL        ", &
                               "N9PMZL        ","N9RAXL        ","N9RAYL        ","N9RAZL        ","N9RDXR        ", &
                               "N9RDYR        ","N9RDZR        ","N9RVXG        ","N9RVYG        ","N9RVZG        ", &
                               "N9TAXL        ","N9TAYL        ","N9TAZL        ","N9TDXR        ","N9TDYR        ", &
                               "N9TDZR        ","N9TVXG        ","N9TVYG        ","N9TVZG        ","PACCACT       ", &
                               "PANGACT       ","PANGINP       ","PRATACT       ","ROOTAPPLIEDFXG","ROOTAPPLIEDFXR", &
                               "ROOTAPPLIEDFYG","ROOTAPPLIEDFYR","ROOTAPPLIEDFZG","ROOTAPPLIEDFZR","ROOTAPPLIEDMXG", &
                               "ROOTAPPLIEDMXR","ROOTAPPLIEDMYG","ROOTAPPLIEDMYR","ROOTAPPLIEDMZG","ROOTAPPLIEDMZR", &
                               "ROOTFXR       ","ROOTFYR       ","ROOTFZR       ","ROOTMXR       ","ROOTMYR       ", &
                               "ROOTMZR       ","TIPRAXL       ","TIPRAYL       ","TIPRAZL       ","TIPRDXR       ", &
                               "TIPRDYR       ","TIPRDZR       ","TIPRVXG       ","TIPRVYG       ","TIPRVZG       ", &
                               "TIPTAXL       ","TIPTAYL       ","TIPTAZL       ","TIPTDXR       ","TIPTDYR       ", &
                               "TIPTDZR       ","TIPTVXG       ","TIPTVYG       ","TIPTVZG       "/)
   INTEGER(IntKi), PARAMETER :: ParamIndxAry(364) =  (/ &                            ! This lists the index into AllOuts(:) of the allowed parameters ValidParamAry(:)
                                        N1DFxl ,         N1DFyl ,         N1DFzl ,         N1DMxl ,         N1DMyl , &
                                        N1DMzl ,          N1Fxl ,          N1Fyl ,          N1Fzl ,          N1Mxl , &
                                         N1Myl ,          N1Mzl ,         N1PFxl ,         N1PFyl ,         N1PFzl , &
                                        N1PMxl ,         N1PMyl ,         N1PMzl ,         N1RAXl ,         N1RAYl , &
                                        N1RAZl ,         N1RDxr ,         N1RDyr ,         N1RDzr ,         N1RVXg , &
                                        N1RVYg ,         N1RVZg ,         N1TAXl ,         N1TAYl ,         N1TAZl , &
                                        N1TDxr ,         N1TDyr ,         N1TDzr ,         N1TVXg ,         N1TVYg , &
                                        N1TVZg ,         N2DFxl ,         N2DFyl ,         N2DFzl ,         N2DMxl , &
                                        N2DMyl ,         N2DMzl ,          N2Fxl ,          N2Fyl ,          N2Fzl , &
                                         N2Mxl ,          N2Myl ,          N2Mzl ,         N2PFxl ,         N2PFyl , &
                                        N2PFzl ,         N2PMxl ,         N2PMyl ,         N2PMzl ,         N2RAXl , &
                                        N2RAYl ,         N2RAZl ,         N2RDxr ,         N2RDyr ,         N2RDzr , &
                                        N2RVXg ,         N2RVYg ,         N2RVZg ,         N2TAXl ,         N2TAYl , &
                                        N2TAZl ,         N2TDxr ,         N2TDyr ,         N2TDzr ,         N2TVXg , &
                                        N2TVYg ,         N2TVZg ,         N3DFxl ,         N3DFyl ,         N3DFzl , &
                                        N3DMxl ,         N3DMyl ,         N3DMzl ,          N3Fxl ,          N3Fyl , &
                                         N3Fzl ,          N3Mxl ,          N3Myl ,          N3Mzl ,         N3PFxl , &
                                        N3PFyl ,         N3PFzl ,         N3PMxl ,         N3PMyl ,         N3PMzl , &
                                        N3RAXl ,         N3RAYl ,         N3RAZl ,         N3RDxr ,         N3RDyr , &
                                        N3RDzr ,         N3RVXg ,         N3RVYg ,         N3RVZg ,         N3TAXl , &
                                        N3TAYl ,         N3TAZl ,         N3TDxr ,         N3TDyr ,         N3TDzr , &
                                        N3TVXg ,         N3TVYg ,         N3TVZg ,         N4DFxl ,         N4DFyl , &
                                        N4DFzl ,         N4DMxl ,         N4DMyl ,         N4DMzl ,          N4Fxl , &
                                         N4Fyl ,          N4Fzl ,          N4Mxl ,          N4Myl ,          N4Mzl , &
                                        N4PFxl ,         N4PFyl ,         N4PFzl ,         N4PMxl ,         N4PMyl , &
                                        N4PMzl ,         N4RAXl ,         N4RAYl ,         N4RAZl ,         N4RDxr , &
                                        N4RDyr ,         N4RDzr ,         N4RVXg ,         N4RVYg ,         N4RVZg , &
                                        N4TAXl ,         N4TAYl ,         N4TAZl ,         N4TDxr ,         N4TDyr , &
                                        N4TDzr ,         N4TVXg ,         N4TVYg ,         N4TVZg ,         N5DFxl , &
                                        N5DFyl ,         N5DFzl ,         N5DMxl ,         N5DMyl ,         N5DMzl , &
                                         N5Fxl ,          N5Fyl ,          N5Fzl ,          N5Mxl ,          N5Myl , &
                                         N5Mzl ,         N5PFxl ,         N5PFyl ,         N5PFzl ,         N5PMxl , &
                                        N5PMyl ,         N5PMzl ,         N5RAXl ,         N5RAYl ,         N5RAZl , &
                                        N5RDxr ,         N5RDyr ,         N5RDzr ,         N5RVXg ,         N5RVYg , &
                                        N5RVZg ,         N5TAXl ,         N5TAYl ,         N5TAZl ,         N5TDxr , &
                                        N5TDyr ,         N5TDzr ,         N5TVXg ,         N5TVYg ,         N5TVZg , &
                                        N6DFxl ,         N6DFyl ,         N6DFzl ,         N6DMxl ,         N6DMyl , &
                                        N6DMzl ,          N6Fxl ,          N6Fyl ,          N6Fzl ,          N6Mxl , &
                                         N6Myl ,          N6Mzl ,         N6PFxl ,         N6PFyl ,         N6PFzl , &
                                        N6PMxl ,         N6PMyl ,         N6PMzl ,         N6RAXl ,         N6RAYl , &
                                        N6RAZl ,         N6RDxr ,         N6RDyr ,         N6RDzr ,         N6RVXg , &
                                        N6RVYg ,         N6RVZg ,         N6TAXl ,         N6TAYl ,         N6TAZl , &
                                        N6TDxr ,         N6TDyr ,         N6TDzr ,         N6TVXg ,         N6TVYg , &
                                        N6TVZg ,         N7DFxl ,         N7DFyl ,         N7DFzl ,         N7DMxl , &
                                        N7DMyl ,         N7DMzl ,          N7Fxl ,          N7Fyl ,          N7Fzl , &
                                         N7Mxl ,          N7Myl ,          N7Mzl ,         N7PFxl ,         N7PFyl , &
                                        N7PFzl ,         N7PMxl ,         N7PMyl ,         N7PMzl ,         N7RAXl , &
                                        N7RAYl ,         N7RAZl ,         N7RDxr ,         N7RDyr ,         N7RDzr , &
                                        N7RVXg ,         N7RVYg ,         N7RVZg ,         N7TAXl ,         N7TAYl , &
                                        N7TAZl ,         N7TDxr ,         N7TDyr ,         N7TDzr ,         N7TVXg , &
                                        N7TVYg ,         N7TVZg ,         N8DFxl ,         N8DFyl ,         N8DFzl , &
                                        N8DMxl ,         N8DMyl ,         N8DMzl ,          N8Fxl ,          N8Fyl , &
                                         N8Fzl ,          N8Mxl ,          N8Myl ,          N8Mzl ,         N8PFxl , &
                                        N8PFyl ,         N8PFzl ,         N8PMxl ,         N8PMyl ,         N8PMzl , &
                                        N8RAXl ,         N8RAYl ,         N8RAZl ,         N8RDxr ,         N8RDyr , &
                                        N8RDzr ,         N8RVXg ,         N8RVYg ,         N8RVZg ,         N8TAXl , &
                                        N8TAYl ,         N8TAZl ,         N8TDxr ,         N8TDyr ,         N8TDzr , &
                                        N8TVXg ,         N8TVYg ,         N8TVZg ,         N9DFxl ,         N9DFyl , &
                                        N9DFzl ,         N9DMxl ,         N9DMyl ,         N9DMzl ,          N9Fxl , &
                                         N9Fyl ,          N9Fzl ,          N9Mxl ,          N9Myl ,          N9Mzl , &
                                        N9PFxl ,         N9PFyl ,         N9PFzl ,         N9PMxl ,         N9PMyl , &
                                        N9PMzl ,         N9RAXl ,         N9RAYl ,         N9RAZl ,         N9RDxr , &
                                        N9RDyr ,         N9RDzr ,         N9RVXg ,         N9RVYg ,         N9RVZg , &
                                        N9TAXl ,         N9TAYl ,         N9TAZl ,         N9TDxr ,         N9TDyr , &
                                        N9TDzr ,         N9TVXg ,         N9TVYg ,         N9TVZg ,        PAccAct , &
                                       PAngAct ,        PAngInp ,        PRatAct , RootAppliedFxg , RootAppliedFxr , &
                                RootAppliedFyg , RootAppliedFyr , RootAppliedFzg , RootAppliedFzr , RootAppliedMxg , &
                                RootAppliedMxr , RootAppliedMyg , RootAppliedMyr , RootAppliedMzg , RootAppliedMzr , &
                                       RootFxr ,        RootFyr ,        RootFzr ,        RootMxr ,        RootMyr , &
                                       RootMzr ,        TipRAXl ,        TipRAYl ,        TipRAZl ,        TipRDxr , &
                                       TipRDyr ,        TipRDzr ,        TipRVXg ,        TipRVYg ,        TipRVZg , &
                                       TipTAXl ,        TipTAYl ,        TipTAZl ,        TipTDxr ,        TipTDyr , &
                                       TipTDzr ,        TipTVXg ,        TipTVYg ,        TipTVZg /)
   CHARACTER(ChanLen), PARAMETER :: ParamUnitsAry(364) =  (/  &  ! This lists the units corresponding to the allowed parameters
                               "(N/m)    ","(N/m)    ","(N/m)    ","(N-m/m)  ","(N-m/m)  ", &
                               "(N-m/m)  ","(N)      ","(N)      ","(N)      ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N)      ","(N)      ","(N)      ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(deg/s^2)","(deg/s^2)", &
                               "(deg/s^2)","(-)      ","(-)      ","(-)      ","(deg/s)  ", &
                               "(deg/s)  ","(deg/s)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ", &
                               "(m)      ","(m)      ","(m)      ","(m/s)    ","(m/s)    ", &
                               "(m/s)    ","(N/m)    ","(N/m)    ","(N/m)    ","(N-m/m)  ", &
                               "(N-m/m)  ","(N-m/m)  ","(N)      ","(N)      ","(N)      ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N)      ","(N)      ", &
                               "(N)      ","(N-m)    ","(N-m)    ","(N-m)    ","(deg/s^2)", &
                               "(deg/s^2)","(deg/s^2)","(-)      ","(-)      ","(-)      ", &
                               "(deg/s)  ","(deg/s)  ","(deg/s)  ","(m/s^2)  ","(m/s^2)  ", &
                               "(m/s^2)  ","(m)      ","(m)      ","(m)      ","(m/s)    ", &
                               "(m/s)    ","(m/s)    ","(N/m)    ","(N/m)    ","(N/m)    ", &
                               "(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N)      ","(N)      ", &
                               "(N)      ","(N-m)    ","(N-m)    ","(N-m)    ","(N)      ", &
                               "(N)      ","(N)      ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ","(-)      ", &
                               "(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(m/s^2)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ","(m)      ", &
                               "(m/s)    ","(m/s)    ","(m/s)    ","(N/m)    ","(N/m)    ", &
                               "(N/m)    ","(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N)      ", &
                               "(N)      ","(N)      ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ", &
                               "(-)      ","(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ", &
                               "(m)      ","(m/s)    ","(m/s)    ","(m/s)    ","(N/m)    ", &
                               "(N/m)    ","(N/m)    ","(N-m/m)  ","(N-m/m)  ","(N-m/m)  ", &
                               "(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(N)      ","(N)      ","(N)      ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(deg/s^2)","(deg/s^2)","(deg/s^2)", &
                               "(-)      ","(-)      ","(-)      ","(deg/s)  ","(deg/s)  ", &
                               "(deg/s)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m)      ", &
                               "(m)      ","(m)      ","(m/s)    ","(m/s)    ","(m/s)    ", &
                               "(N/m)    ","(N/m)    ","(N/m)    ","(N-m/m)  ","(N-m/m)  ", &
                               "(N-m/m)  ","(N)      ","(N)      ","(N)      ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N)      ","(N)      ","(N)      ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(deg/s^2)","(deg/s^2)", &
                               "(deg/s^2)","(-)      ","(-)      ","(-)      ","(deg/s)  ", &
                               "(deg/s)  ","(deg/s)  ","(m/s^2)  ","(m/s^2)  ","(m/s^2)  ", &
                               "(m)      ","(m)      ","(m)      ","(m/s)    ","(m/s)    ", &
                               "(m/s)    ","(N/m)    ","(N/m)    ","(N/m)    ","(N-m/m)  ", &
                               "(N-m/m)  ","(N-m/m)  ","(N)      ","(N)      ","(N)      ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N)      ","(N)      ", &
                               "(N)      ","(N-m)    ","(N-m)    ","(N-m)    ","(deg/s^2)", &
                               "(deg/s^2)","(deg/s^2)","(-)      ","(-)      ","(-)      ", &
                               "(deg/s)  ","(deg/s)  ","(deg/s)  ","(m/s^2)  ","(m/s^2)  ", &
                               "(m/s^2)  ","(m)      ","(m)      ","(m)      ","(m/s)    ", &
                               "(m/s)    ","(m/s)    ","(N/m)    ","(N/m)    ","(N/m)    ", &
                               "(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N)      ","(N)      ", &
                               "(N)      ","(N-m)    ","(N-m)    ","(N-m)    ","(N)      ", &
                               "(N)      ","(N)      ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ","(-)      ", &
                               "(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ","(m/s^2)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ","(m)      ", &
                               "(m/s)    ","(m/s)    ","(m/s)    ","(N/m)    ","(N/m)    ", &
                               "(N/m)    ","(N-m/m)  ","(N-m/m)  ","(N-m/m)  ","(N)      ", &
                               "(N)      ","(N)      ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ", &
                               "(-)      ","(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ", &
                               "(m)      ","(m/s)    ","(m/s)    ","(m/s)    ","(deg/s^2)", &
                               "(deg)    ","(deg)    ","(deg/s)  ","(N)      ","(N)      ", &
                               "(N)      ","(N)      ","(N)      ","(N)      ","(N-m)    ", &
                               "(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ","(N-m)    ", &
                               "(N)      ","(N)      ","(N)      ","(N-m)    ","(N-m)    ", &
                               "(N-m)    ","(deg/s^2)","(deg/s^2)","(deg/s^2)","(-)      ", &
                               "(-)      ","(-)      ","(deg/s)  ","(deg/s)  ","(deg/s)  ", &
                               "(m/s^2)  ","(m/s^2)  ","(m/s^2)  ","(m)      ","(m)      ", &
                               "(m)      ","(m/s)    ","(m/s)    ","(m/s)    "/)

   INTEGER, PARAMETER :: RootAppliedLd(12) = (/ &  ! all applied load mapped to root outputs.  For convenience in setting calculation flag
                              RootAppliedFxr, RootAppliedFyr, RootAppliedFzr, RootAppliedMxr, RootAppliedMyr, RootAppliedMzr, &
                              RootAppliedFxg, RootAppliedFyg, RootAppliedFzg, RootAppliedMxg, RootAppliedMyg, RootAppliedMzg /)


      ! Initialize values
   ErrStat = ErrID_None
   ErrMsg = ""
   InvalidOutput = .FALSE.


!   ..... Developer must add checking for invalid inputs here: .....
   DO I = p%NNodeOuts+1,9  ! Invalid nodes

      InvalidOutput( NFl( i,:) ) = .true.
      InvalidOutput( NMl( i,:) ) = .true.
      InvalidOutput( NTDr(i,:) ) = .true.
      InvalidOutput( NRDr(i,:) ) = .true.
      InvalidOutput( NTVg(i,:) ) = .true.
      InvalidOutput( NRVg(i,:) ) = .true.
      InvalidOutput( NTAl(i,:) ) = .true.
      InvalidOutput( NRAl(i,:) ) = .true.
      InvalidOutput( NPFl(i,:) ) = .true.
      InvalidOutput( NPMl(i,:) ) = .true.
      InvalidOutput( NDFl(i,:) ) = .true.
      InvalidOutput( NDMl(i,:) ) = .true.
         
   END DO   
   
   IF (.not. p%UsePitchAct) THEN
      InvalidOutput( PAngInp ) = .true.
      InvalidOutput( PAngAct ) = .true.
      InvalidOutput( PRatAct ) = .true.
      InvalidOutput( PAccAct ) = .true.
   END IF

   if (p%BldMotionNodeLoc /= BD_MESH_FE) then
      DO I = 1,9
         InvalidOutput( NPFl(i,:) ) = .true.
         InvalidOutput( NPMl(i,:) ) = .true. 
      END DO      
   end if

!   ................. End of validity checking .................


   !-------------------------------------------------------------------------------------------------
   ! Allocate and set index, name, and units for the output channels
   ! If a selected output channel is not available in this module, set error flag.
   !-------------------------------------------------------------------------------------------------

   ALLOCATE ( p%OutParam(0:p%NumOuts) , STAT=ErrStat2 )
   IF ( ErrStat2 /= 0_IntKi )  THEN
      CALL SetErrStat( ErrID_Fatal,"Error allocating memory for the BeamDyn OutParam array.", ErrStat, ErrMsg, RoutineName )
      RETURN
   ENDIF

      ! Set index, name, and units for the time output channel:

   p%OutParam(0)%Indx  = Time
   p%OutParam(0)%Name  = "Time"    ! OutParam(0) is the time channel by default.
   p%OutParam(0)%Units = "(s)"
   p%OutParam(0)%SignM = 1


      ! Set index, name, and units for all of the output channels.
      ! If a selected output channel is not available by this module set ErrStat = ErrID_Warn.

   DO I = 1,p%NumOuts

      p%OutParam(I)%Name  = OutList(I)
      OutListTmp          = OutList(I)

      ! Reverse the sign (+/-) of the output channel if the user prefixed the
      !   channel name with a "-", "_", "m", or "M" character indicating "minus".


      CheckOutListAgain = .FALSE.

      IF      ( INDEX( "-_", OutListTmp(1:1) ) > 0 ) THEN
         p%OutParam(I)%SignM = -1                         ! ex, "-TipDxc1" causes the sign of TipDxc1 to be switched.
         OutListTmp          = OutListTmp(2:)
      ELSE IF ( INDEX( "mM", OutListTmp(1:1) ) > 0 ) THEN ! We'll assume this is a variable name for now, (if not, we will check later if OutListTmp(2:) is also a variable name)
         CheckOutListAgain   = .TRUE.
         p%OutParam(I)%SignM = 1
      ELSE
         p%OutParam(I)%SignM = 1
      END IF

      CALL Conv2UC( OutListTmp )    ! Convert OutListTmp to upper case


      Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )


         ! If it started with an "M" (CheckOutListAgain) we didn't find the value in our list (Indx < 1)

      IF ( CheckOutListAgain .AND. Indx < 1 ) THEN    ! Let's assume that "M" really meant "minus" and then test again
         p%OutParam(I)%SignM = -1                     ! ex, "MTipDxc1" causes the sign of TipDxc1 to be switched.
         OutListTmp          = OutListTmp(2:)

         Indx = IndexCharAry( OutListTmp(1:OutStrLenM1), ValidParamAry )
      END IF


      IF ( Indx > 0 ) THEN ! we found the channel name
         p%OutParam(I)%Indx     = ParamIndxAry(Indx)
         IF ( InvalidOutput( ParamIndxAry(Indx) ) ) THEN  ! but, it isn't valid for these settings
            p%OutParam(I)%Units = "INVALID"
            p%OutParam(I)%SignM = 0
         ELSE
            p%OutParam(I)%Units = ParamUnitsAry(Indx) ! it's a valid output
            
            if ( p%OutParam(I)%Indx >= N1DFxl .and. p%OutParam(I)%Indx <= N9DMzl ) p%OutInputs = .true.
            if ( any(RootAppliedLd==p%OutParam(I)%Indx) ) p%CompAppliedLdAtRoot = .true.  ! need to setup meshes
         END IF
      ELSE ! this channel isn't valid
         p%OutParam(I)%Indx  = Time                 ! pick any valid channel (I just picked "Time" here because it's universal)
         p%OutParam(I)%Units = "INVALID"
         p%OutParam(I)%SignM = 0                    ! multiply all results by zero

         CALL SetErrStat(ErrID_Fatal, TRIM(p%OutParam(I)%Name)//" is not an available output channel.",ErrStat,ErrMsg,RoutineName)
      END IF

   END DO

   RETURN
END SUBROUTINE SetOutParam
!----------------------------------------------------------------------------------------------------------------------------------
!End of code generated by Matlab script
!**********************************************************************************************************************************

!----------------------------------------------------------------------------------------------------------------------------------
!> This routine validates the inputs from the BeamDyn input files.
SUBROUTINE BD_ValidateInputData( InitInp, InputFileData, ErrStat, ErrMsg )

      ! Passed variables:

   TYPE(BD_InitInputType),  INTENT(IN   )  :: InitInp           !< Input data for initialization routine
   TYPE(BD_InputFile),      INTENT(IN   )  :: InputFileData     !< All the data in the BeamDyn input file
   INTEGER(IntKi),          INTENT(  OUT)  :: ErrStat           !< Error status
   CHARACTER(*),            INTENT(  OUT)  :: ErrMsg            !< Error message


      ! local variables
   INTEGER(IntKi)                          :: i,j                                 ! loop counters
   INTEGER(IntKi)                          :: nNodes                              ! number of nodes that will be on the BldMotion mesh
   REAL(SiKi)                              :: r1, r2                              ! use single-precision real numbers for validity check

   CHARACTER(*), PARAMETER                 :: RoutineName = 'BD_ValidateInputData'

   ErrStat = ErrID_None
   ErrMsg  = ""


      ! QuasiStaticInit only valid with DynamicSolve.  Warn and ignore
   IF(InputFileData%QuasiStaticInit .AND. (.NOT. InitInp%DynamicSolve)) & 
      CALL SetErrStat ( ErrID_Warn, 'QuasiStaticInit option only valid with dynamic solve.  Ignoring.', ErrStat, ErrMsg, RoutineName )

   IF(InputFileData%rhoinf .LT. 0.0_BDKi .OR. InputFileData%rhoinf .GT. 1.0_BDKi) &
       CALL SetErrStat ( ErrID_Fatal, 'Numerical damping parameter \rho_{inf} must be in the range of [0.0,1.0]', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%quadrature .NE. GAUSS_QUADRATURE .AND. InputFileData%quadrature .NE. TRAP_QUADRATURE) &
       CALL SetErrStat ( ErrID_Fatal, 'Quadrature type must be 1 (Gauss) or 2 (Trapezoidal)', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%refine .LT. 1 ) &
       CALL SetErrStat ( ErrID_Fatal, 'Refinement parameter must be greater than or equal to 1', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%n_fact .LT. 1 ) &
       CALL SetErrStat ( ErrID_Fatal, 'Factorization parameter must be greater than or equal to 1', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%member_total .LT. 1 ) &
       CALL SetErrStat ( ErrID_Fatal, 'member_total must be greater than 0', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%member_total .NE. 1 .AND. InputFileData%quadrature .EQ. TRAP_QUADRATURE) &
       CALL SetErrStat ( ErrID_Fatal, 'Trapzoidal quadrature only allows one member (element)', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%kp_total .LT. 3 ) THEN
       CALL SetErrStat ( ErrID_Fatal, 'kp_total must be greater than or equal to 3', ErrStat, ErrMsg, RoutineName )
   ELSE ! added this in the "else" in case InputFileData%kp_coordinate isn't allocated with at least 1 point
      IF ( .not. EqualRealNos( InputFileData%kp_coordinate(1,3), 0.0_BDKi ) ) THEN
         CALL SetErrStat(ErrID_Fatal, 'kp_zr on first key point must be 0.', ErrStat, ErrMsg, RoutineName )
      END IF
         !bjj: added checks on kp_xr(1) and kp_yr(1) because BeamDyn's equations currently assume that the blade root and the first FE node 
         !     are at the same point (i.e., u%BladeRoot is located at the first finite-element node)
      IF ( .not. EqualRealNos( InputFileData%kp_coordinate(1,2), 0.0_BDKi ) ) THEN
         CALL SetErrStat(ErrID_Fatal, 'kp_yr on first key point must be 0.', ErrStat, ErrMsg, RoutineName )
      END IF
      IF ( .not. EqualRealNos( InputFileData%kp_coordinate(1,1), 0.0_BDKi ) ) THEN
         CALL SetErrStat(ErrID_Fatal, 'kp_xr on first key point must be 0.', ErrStat, ErrMsg, RoutineName )
      END IF
   END IF

   DO i=1,InputFileData%member_total
       IF(InputFileData%kp_member(i) .LT. 3) THEN
          CALL SetErrStat(ErrID_Fatal,'There must be at least three key points in member number '//TRIM(Num2LStr(i))//'.', ErrStat, ErrMsg,RoutineName)
          EXIT
       ENDIF
   ENDDO
   IF(SUM(InputFileData%kp_member) .NE. InputFileData%kp_total+InputFileData%member_total-1 ) &
       CALL SetErrStat ( ErrID_Fatal, 'Geometric definition error: kp_total and key points in each member are inconsistent', ErrStat, ErrMsg,RoutineName)

   IF(InputFileData%order_elem .LT. 1 ) &
       CALL SetErrStat ( ErrID_Fatal, 'order_elem must be greater than 0', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%order_elem .LT. 2 .AND. InputFileData%quadrature .EQ. GAUSS_QUADRATURE) &
       CALL SetErrStat ( ErrID_Fatal, 'order_elem must be greater than 1 for Gauss quadrature', ErrStat, ErrMsg, RoutineName )

   IF(InputFileData%InpBl%station_total .LT. 2 ) &
       CALL SetErrStat ( ErrID_Fatal, 'Number of material stations along blade axis much be greater than 2', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%InpBl%station_eta(1) .NE. 0.0_BDKi ) &
       CALL SetErrStat ( ErrID_Fatal, 'The first station_eta must be equal to 0.0 (root)', ErrStat, ErrMsg, RoutineName )
   IF(InputFileData%InpBl%station_eta(InputFileData%InpBl%station_total) .NE. 1.0_BDKi ) &
       CALL SetErrStat ( ErrID_Fatal, 'The last station_eta must be equal to 1.0 (tip)', ErrStat, ErrMsg, RoutineName )

   IF (InputFileData%NRMax < 2) CALL SetErrStat ( ErrID_Fatal, 'Maximum number of iterations in Newton-Raphson (NRMax) must be greater than 1.', ErrStat, ErrMsg, RoutineName )
   IF (InputFileData%stop_tol < EPSILON(InputFileData%stop_tol) ) &
      CALL SetErrStat ( ErrID_Fatal, 'Tolerance for stopping (stop_tol) must be larger than machine precision ('//trim(num2lstr(EPSILON(InputFileData%stop_tol)))//').', ErrStat, ErrMsg, RoutineName )

   if (InputFileData%UsePitchAct) then
      if ( EqualRealNos(InputFileData%pitchJ, 0.0_BDKi) ) call SetErrStat(ErrID_Fatal,'Pitch actuator inertia must not be 0.',ErrStat,ErrMsg,RoutineName)
   end if

   !bjj: should check that there aren't any overlapping trapezoidal quadrature points (is this done in the input file section??)

      ! Check that the values for the damping are similar.
      !  According to Qi, the damping values mu1 .. mu6 should be of the same order of magnitude.  If they aren't, BeamDyn will likely fail badly.
      !  Will assume for moment that they must be within a factor of 5 of the first value given.
!FIXME: Check a valid range sometime.
   IF ( (maxval(abs(InputFileData%InpBl%beta)) / minval(abs(InputFileData%InpBl%beta)) > 10.0_R8Ki) .AND. InputFileData%InpBl%damp_flag == 1 ) THEN
      call SetErrStat( ErrID_Warn,'Damping values in blade file are not of similar order of magnitude. BeamDyn may not converge!', ErrStat, ErrMsg, RoutineName )
   ENDIF


   DO j=1,InputFileData%InpBl%station_total
      DO i=2,3
         IF ( .not. EqualRealNos( InputFileData%InpBl%mass0(i,i,j), InputFileData%InpBl%mass0(1,1,j) ) ) then
            call SetErrStat( ErrID_Fatal, 'Input station '//trim(num2lstr(j))//' mass densities are not the same (i.e., first 3 diagonal elements in mass matrix must be equal).', ErrStat, ErrMsg, RoutineName )
            exit
         END IF
      END DO

      !............
      ! NOTE: InputFileData%InpBl%mass0 is in IEC coordinates (it was previously in beam coordinates); error message refers to IEC coordinates in input file
      r1 = InputFileData%InpBl%mass0(6,6,j)
      r2 = InputFileData%InpBl%mass0(4,4,j) + InputFileData%InpBl%mass0(5,5,j)
      IF ( abs(r1 - r2) > Mass0Tol ) then
         call SetErrStat( ErrID_Fatal, 'Input station '//trim(num2lstr(j))//' i_plr must equal i_Edg + i_Flp (i.e., sum of 4th and 5th diagonal elements in mass matrix must equal the 6th diagonal element).', ErrStat, ErrMsg, RoutineName )
      END IF
      !............
   END DO


      ! .............................
      ! check outputs:
      ! .............................

   if ( ( InputFileData%NNodeOuts < 0_IntKi ) .OR. ( InputFileData%NNodeOuts > 9_IntKi ) )  then
      call SetErrStat( ErrID_Fatal, 'NNodeOuts must be between 0 and 9 (inclusive).', ErrStat, ErrMsg, RoutineName )
   else

   ! Check to see if all OutNd(:) analysis points are existing analysis points:
!bjj: FIXME!!!! THIS ISN'T NECESSARIALLY THE NUMBER OF NODES ON THIS MESH
!      IF (p%BldMotionNodeLoc == BD_MESH_QP) THEN
         if(InputFileData%quadrature .EQ. TRAP_QUADRATURE) then
            nNodes = (InputFileData%InpBl%station_total - 1)*InputFileData%refine + 1  ! number of nodes on y%BldMotion mesh
         else
            nNodes = (InputFileData%order_elem + 1)*InputFileData%member_total  ! = p%nodes_per_elem*p%elem_total (number of nodes on y%BldMotion mesh)
         end if
!      ELSE
!         nNodes = (InputFileData%order_elem + 1)*InputFileData%member_total + 1
!      END IF

      do j=1,InputFileData%NNodeOuts
         if ( InputFileData%OutNd(j) < 1_IntKi .OR. InputFileData%OutNd(j) > nNodes ) then
            call SetErrStat( ErrID_Fatal, ' All OutNd values must be between 1 and '//&
                    trim( Num2LStr( nNodes ) )//' (inclusive).', ErrStat, ErrMsg, RoutineName )
            exit ! stop checking this loop
         end if
      end do

   end if

   !..................
   ! check for linearization
   !..................
   if (InitInp%Linearize) then
      if (.NOT. InitInp%DynamicSolve) THEN
         call SetErrStat( ErrID_Fatal, 'Static analysis cannot be used for BeamDyn linearization. Set DynamicSolve to TRUE.', ErrStat, ErrMsg, RoutineName )
      end if
      
      if (InputFileData%UsePitchAct) then
         call SetErrStat( ErrID_Fatal, 'Pitch actuator model cannot currently be used for linearization in BeamDyn. Set UsePitchAct=False.', ErrStat, ErrMsg, RoutineName )
      end if
   end if
   
END SUBROUTINE BD_ValidateInputData
!----------------------------------------------------------------------------------------------------------------------------------
!> this routine fills the AllOuts array, which is used to send data to the glue code to be written to an output file.
SUBROUTINE Calc_WriteOutput( p, AllOuts, y, m, ErrStat, ErrMsg, CalcWriteOutput )

   TYPE(BD_ParameterType),       INTENT(IN   )  :: p                                 !< The module parameters
   REAL(ReKi),                   INTENT(INOUT)  :: AllOuts(0:)                       !< array of values to potentially write to file
   TYPE(BD_OutputType),          INTENT(IN   )  :: y                                 !< outputs
   TYPE(BD_MiscVarType),         INTENT(INOUT)  :: m                                 !< misc/optimization variables (for computing mesh transfers)
   INTEGER(IntKi),               INTENT(  OUT)  :: ErrStat                           !< The error status code
   CHARACTER(*),                 INTENT(  OUT)  :: ErrMsg                            !< The error message, if an error occurred
   LOGICAL               ,       INTENT(IN   )  :: CalcWriteOutput                   !< flag that determines if we need to compute AllOuts (or just the reaction loads that get returned to ServoDyn)

      ! local variables
   CHARACTER(*), PARAMETER                      :: RoutineName = 'Calc_WriteOutput'
   INTEGER(IntKi)                               :: ErrStat2
   CHARACTER(ErrMsgLen)                         :: ErrMsg2


   INTEGER(IntKi)                               :: j,beta,j_BldMotion
   REAL(BDKi)                                   :: temp_vec(3)
   REAL(BDKi)                                   :: temp_vec2(3)
   REAL(BDKi)                                   :: temp33(3,3)
   REAL(BDKi)                                   :: temp33_2(3,3)
   
   REAL(BDKi)                                   :: d_ref(3)
   REAL(BDKi)                                   :: d(3)
   REAL(BDKi)                                   :: RootRelOrient(3,3)

   

      ! start routine:
   ErrStat = ErrID_None
   ErrMsg  = ""

      !-------------------------
      ! Reaction forces
   temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%ReactionForce%Force(:,1))
   AllOuts( RootFxr ) = temp_vec(1)
   AllOuts( RootFyr ) = temp_vec(2)
   AllOuts( RootFzr ) = temp_vec(3)

      !-------------------------
      ! Reaction moments (these are being computed for ServoDyn output, as well as for WriteOutput)
   temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),y%ReactionForce%Moment(:,1))
   AllOuts( RootMxr ) = temp_vec(1)
   AllOuts( RootMyr ) = temp_vec(2)
   AllOuts( RootMzr ) = temp_vec(3)

   !-------------------------
   ! we don't need to calculate the rest of these values if we don't ask for WriteOutput channels
   ! (but we did need RootMxr and RootMyr)
   if ( p%NumOuts <= 0 .or. .not. CalcWriteOutput) RETURN
   !-------------------------

   
      ! compute the root relative orientation, RootRelOrient, which is used in several calculations below 
      ! RootRelOrient = matmul( transpose(m%u2%RootMotion%Orientation(:,:,1)), m%u2%RootMotion%RefOrientation(:,:,1))
   call LAPACK_GEMM('T', 'N', 1.0_BDKi, m%u2%RootMotion%Orientation(:,:,1), m%u2%RootMotion%RefOrientation(:,:,1), 0.0_BDKi, RootRelOrient,   ErrStat2, ErrMsg2 )
   
      !------------------------------------
      ! Tip translational deflection (relative to the undeflected position) expressed in r   
   d     = y%BldMotion%TranslationDisp(:, y%BldMotion%NNodes) - m%u2%RootMotion%TranslationDisp(:,1)
   d_ref = y%BldMotion%Position(       :, y%BldMotion%NNodes) - m%u2%RootMotion%Position(       :,1)
   temp_vec2 = d + d_ref - matmul( RootRelOrient, d_ref ) ! tip displacement
   temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2)
      
   AllOuts( TipTDxr ) = temp_vec(1)
   AllOuts( TipTDyr ) = temp_vec(2)
   AllOuts( TipTDzr ) = temp_vec(3)

      !-------------------------
      ! Tip angular/rotational deflection Wiener-Milenkovic parameter (relative to the undeflected orientation) expressed in r
   call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,y%BldMotion%NNodes), RootRelOrient,   0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 )
   call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation(   :,:,y%BldMotion%NNodes), temp33_2,        0.0_BDKi, temp33,   ErrStat2, ErrMsg2 )
   call BD_CrvExtractCrv(temp33,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the tip angular/rotational defelctions
      CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
      if (ErrStat >= AbortErrLev) return
   temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2) ! translate these parameters to the correct system for output
   
   AllOuts( TipRDxr ) = temp_vec(1)
   AllOuts( TipRDyr ) = temp_vec(2)
   AllOuts( TipRDzr ) = temp_vec(3)

      !-------------------------
      ! Tip translational velocities (absolute) expressed in g, given in m/s
   AllOuts( TipTVXg ) = y%BldMotion%TranslationVel(1,y%BldMotion%NNodes)
   AllOuts( TipTVYg ) = y%BldMotion%TranslationVel(2,y%BldMotion%NNodes)
   AllOuts( TipTVZg ) = y%BldMotion%TranslationVel(3,y%BldMotion%NNodes)

      !-------------------------
      ! Tip  angular/rotational velocities (absolute) expressed in g, given in deg/s
   AllOuts( TipRVXg ) = y%BldMotion%RotationVel(1,y%BldMotion%NNodes)*R2D
   AllOuts( TipRVYg ) = y%BldMotion%RotationVel(2,y%BldMotion%NNodes)*R2D
   AllOuts( TipRVZg ) = y%BldMotion%RotationVel(3,y%BldMotion%NNodes)*R2D

      !-------------------------
      ! Tip translational accelerations (absolute) expressed in l, given in m/s^2
   temp_vec = MATMUL(y%BldMotion%Orientation(:,:,y%BldMotion%NNodes), y%BldMotion%TranslationAcc(:,y%BldMotion%NNodes)) ! translate accelerations to local system for output   
   AllOuts( TipTAXl ) = temp_vec(1)
   AllOuts( TipTAYl ) = temp_vec(2)
   AllOuts( TipTAZl ) = temp_vec(3)

      !-------------------------
      ! Tip angular/rotational accelerations (absolute) expressed in l, given in deg/s^2
   temp_vec = MATMUL(y%BldMotion%Orientation(:,:,y%BldMotion%NNodes), y%BldMotion%RotationAcc(:,y%BldMotion%NNodes)) ! translate accelerations to local system for output   
   AllOuts( TipRAXl ) = temp_vec(1)*R2D
   AllOuts( TipRAYl ) = temp_vec(2)*R2D
   AllOuts( TipRAZl ) = temp_vec(3)*R2D


   !------------------------------------
   ! Calculate Motion of blade nodes
   !------------------------------------

      ! outputs on the nodes
   DO beta=1,p%NNodeOuts

      j=p%OutNd(beta)
      j_BldMotion = p%NdIndx(j)

         !------------------------------------
         ! Sectional force resultants at Node 1 expressed in l, given in N
      SELECT CASE (p%BldMotionNodeLoc)
      CASE (BD_MESH_FE)
         temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceFE(1:3,j))
      CASE (BD_MESH_QP)
         temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceQP(1:3,j))
      END SELECT
      AllOuts( NFl( beta,1 ) ) = temp_vec(1)
      AllOuts( NFl( beta,2 ) ) = temp_vec(2)
      AllOuts( NFl( beta,3 ) ) = temp_vec(3)
      
         !------------------------------------
         ! Sectional moment resultants at Node 1 expressed in l, given in N-m
      SELECT CASE (p%BldMotionNodeLoc)
      CASE (BD_MESH_FE)
         temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceFE(4:6,j))
      CASE (BD_MESH_QP)
         temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%BldInternalForceQP(4:6,j))
      END SELECT
      AllOuts( NMl( beta,1 ) ) = temp_vec(1)
      AllOuts( NMl( beta,2 ) ) = temp_vec(2)
      AllOuts( NMl( beta,3 ) ) = temp_vec(3)
         
      
         !------------------------------------
         ! Sectional translational deflection (relative to the undeflected position) expressed in r   
      d     = y%BldMotion%TranslationDisp(:, j_BldMotion) - m%u2%RootMotion%TranslationDisp(:,1)
      d_ref = y%BldMotion%Position(       :, j_BldMotion) - m%u2%RootMotion%Position(       :,1)
      temp_vec2 = d + d_ref - MATMUL( RootRelOrient, d_ref ) ! tip displacement
      temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2)
      
      AllOuts( NTDr( beta,1 ) ) = temp_vec(1)
      AllOuts( NTDr( beta,2 ) ) = temp_vec(2)
      AllOuts( NTDr( beta,3 ) ) = temp_vec(3)

         !-------------------------
         ! Sectional angular/rotational deflection Wiener-Milenkovic parameter (relative to the undeflected orientation) expressed in r
      call LAPACK_GEMM('N', 'T', 1.0_BDKi, y%BldMotion%RefOrientation(:,:,j_BldMotion), RootRelOrient,  0.0_BDKi, temp33_2, ErrStat2, ErrMsg2 )
      call LAPACK_GEMM('T', 'N', 1.0_BDKi, y%BldMotion%Orientation(   :,:,j_BldMotion), temp33_2,       0.0_BDKi, temp33,   ErrStat2, ErrMsg2 )
      call BD_CrvExtractCrv(temp33,temp_vec2, ErrStat2, ErrMsg2) ! temp_vec2 = the Wiener-Milenkovic parameters of the node's angular/rotational defelctions
         CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
         if (ErrStat >= AbortErrLev) return
      temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),temp_vec2) ! translate these parameters to the correct system for output
            
      AllOuts( NRDr( beta,1 ) ) = temp_vec(1)
      AllOuts( NRDr( beta,2 ) ) = temp_vec(2)
      AllOuts( NRDr( beta,3 ) ) = temp_vec(3)

         !-------------------------
         ! Sectional translational velocities (absolute) expressed in g, given in m/s
      AllOuts( NTVg( beta,1 ) ) = y%BldMotion%TranslationVel(1,j_BldMotion)
      AllOuts( NTVg( beta,2 ) ) = y%BldMotion%TranslationVel(2,j_BldMotion)
      AllOuts( NTVg( beta,3 ) ) = y%BldMotion%TranslationVel(3,j_BldMotion)

         !-------------------------
         ! Sectional  angular/rotational velocities (absolute) expressed in g, given in deg/s
      AllOuts( NRVg( beta,1 ) ) = y%BldMotion%RotationVel(1,j_BldMotion)*R2D
      AllOuts( NRVg( beta,2 ) ) = y%BldMotion%RotationVel(2,j_BldMotion)*R2D
      AllOuts( NRVg( beta,3 ) ) = y%BldMotion%RotationVel(3,j_BldMotion)*R2D


         !-------------------------
         ! Sectional Tip translational accelerations (absolute) expressed in l, given in m/s^2
      temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), y%BldMotion%TranslationAcc(:,j_BldMotion))      
      AllOuts( NTAl( beta,1 ) ) = temp_vec(1)
      AllOuts( NTAl( beta,2 ) ) = temp_vec(2)
      AllOuts( NTAl( beta,3 ) ) = temp_vec(3)

         !-------------------------
         ! Sectional angular/rotational accelerations (absolute) expressed in l, given in deg/s^2
      temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), y%BldMotion%RotationAcc(:,j_BldMotion))      
      AllOuts( NRAl( beta,1 ) ) = temp_vec(1)*R2D
      AllOuts( NRAl( beta,2 ) ) = temp_vec(2)*R2D
      AllOuts( NRAl( beta,3 ) ) = temp_vec(3)*R2D


      if (p%BldMotionNodeLoc == BD_MESH_FE) THEN !> FIXME: If we are on the finite element points, the input and output meshes are siblings, otherwise we need to multiply by a different orientation (if we're okay 
                                                 !! with the nodes meaning something different) or we need to map the u2%PointLoad like we do for the m%u2%DistrLoad%Force loads.
            !-------------------------
            ! Applied point forces at Node 1 expressed in l, given in N
         temp_vec = MATMUL(y%BldMotion%Orientation(1:3,1:3,j_BldMotion),m%u2%PointLoad%Force(:,j))
         AllOuts( NPFl( beta,1 ) ) = temp_vec(1)
         AllOuts( NPFl( beta,2 ) ) = temp_vec(2)
         AllOuts( NPFl( beta,3 ) ) = temp_vec(3)

            !-------------------------
            ! Applied point moments at Node 1 expressed in l, given in N-m
         temp_vec = MATMUL(y%BldMotion%Orientation(1:3,1:3,j_BldMotion),m%u2%PointLoad%Moment(:,j))
         AllOuts( NPMl( beta,1 ) ) = temp_vec(1)
         AllOuts( NPMl( beta,2 ) ) = temp_vec(2)
         AllOuts( NPMl( beta,3 ) ) = temp_vec(3)
      end if
      

   end do ! nodes

      ! to avoid unnecessary mesh mapping calculations, calculate these outputs only when we've requested them
   if (p%OutInputs) then
               
      if (p%BldMotionNodeLoc == BD_MESH_QP) THEN ! If we are on the quadrature points, the input and output meshes are siblings
         
         do beta=1,p%NNodeOuts

            j=p%OutNd(beta)
            j_BldMotion = p%NdIndx(j)

               !-------------------------
               ! Applied distributed forces at Node 1 expressed in l, in N/m
            temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%u2%DistrLoad%Force( :,j_BldMotion))
            AllOuts( NDFl( beta,1 ) ) = temp_vec(1)
            AllOuts( NDFl( beta,2 ) ) = temp_vec(2)
            AllOuts( NDFl( beta,3 ) ) = temp_vec(3)
         
               !-------------------------
               ! Applied distributed moments at Node 1 expressed in l, in N-m/m
            temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%u2%DistrLoad%Moment(:,j_BldMotion))
            AllOuts( NDMl( beta,1 ) ) = temp_vec(1)
            AllOuts( NDMl( beta,2 ) ) = temp_vec(2)
            AllOuts( NDMl( beta,3 ) ) = temp_vec(3)

         end do ! nodes
         
         
      else
         
         
            ! transfer the output motions to the input nodes for load transfer
         CALL Transfer_Line2_to_Line2( y%BldMotion, m%y_BldMotion_at_u, m%Map_y_BldMotion_to_u, ErrStat2, ErrMsg2 )
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

            ! transfer the input loads to the output nodes for writing output
         CALL Transfer_Line2_to_Line2( m%u2%DistrLoad, m%u_DistrLoad_at_y, m%Map_u_DistrLoad_to_y, ErrStat2, ErrMsg2, m%y_BldMotion_at_u, y%BldMotion)
            CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )


         do beta=1,p%NNodeOuts

            j=p%OutNd(beta)
            j_BldMotion = p%NdIndx(j)

               !-------------------------
               ! Applied distributed forces at Node 1 expressed in l, in N/m
            temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%u_DistrLoad_at_y%Force( :,j_BldMotion))
            AllOuts( NDFl( beta,1 ) ) = temp_vec(1)
            AllOuts( NDFl( beta,2 ) ) = temp_vec(2)
            AllOuts( NDFl( beta,3 ) ) = temp_vec(3)
         
               !-------------------------
               ! Applied distributed moments at Node 1 expressed in l, in N-m/m
            temp_vec = MATMUL(y%BldMotion%Orientation(:,:,j_BldMotion), m%u_DistrLoad_at_y%Moment(:,j_BldMotion))
            AllOuts( NDMl( beta,1 ) ) = temp_vec(1)
            AllOuts( NDMl( beta,2 ) ) = temp_vec(2)
            AllOuts( NDMl( beta,3 ) ) = temp_vec(3)

         end do ! nodes
         
      end if
   end if

   ! compute mapping of applied distributed loads to the root location
   !   FIXME: currently not mapping the PointLoads over since there is no motion mesh associated with the PointLoad
   !        To get the PointLoads mapped, the following is necessary
   !        1. create a y%BldMotionFE mesh at the finite element points (this would be very useful someday)
   !           - take position info directly from uuu and state information etc.
   !           - Add output channels for this for comparison
   !        2. make u%PointLoad a sibling of y%BldMotionFE
   !        3. Setup m%Map_u_PtLoad_to_R
   !        4. map loads and do summation here (remember Transfer zero's out forces/moments, so add temp arrays to hold those)
   if (p%CompAppliedLdAtRoot .and. p%BldMotionNodeLoc == BD_MESH_QP) then

      ! shorthand to simplify life
      associate(RF => m%LoadsAtRoot%Force, RM => m%LoadsAtRoot%Moment)

      ! mapping of distributed loads to LoadsAtRoot
      call Transfer_Line2_to_Point( m%u2%DistrLoad, m%LoadsAtRoot, m%Map_u_DistrLoad_to_R, ErrStat2, ErrMsg2, y%BldMotion, m%u2%RootMotion )
         call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

      ! Global coords
      AllOuts( RootAppliedFxg ) = RF(1,1)
      AllOuts( RootAppliedFyg ) = RF(2,1)
      AllOuts( RootAppliedFzg ) = RF(3,1)
      AllOuts( RootAppliedMxg ) = RM(1,1)
      AllOuts( RootAppliedMyg ) = RM(2,1)
      AllOuts( RootAppliedMzg ) = RM(3,1)

      ! Root coords
      temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),RF(:,1))
      AllOuts( RootAppliedFxr ) = temp_vec(1)
      AllOuts( RootAppliedFyr ) = temp_vec(2)
      AllOuts( RootAppliedFzr ) = temp_vec(3)
      temp_vec = MATMUL(m%u2%RootMotion%Orientation(:,:,1),RM(:,1))
      AllOuts( RootAppliedMxr ) = temp_vec(1)
      AllOuts( RootAppliedMyr ) = temp_vec(2)
      AllOuts( RootAppliedMzr ) = temp_vec(3)

      end associate
   endif


END SUBROUTINE Calc_WriteOutput
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine generates the summary file, which contains a regurgitation of  the input data and interpolated flexible body data.
SUBROUTINE BD_PrintSum( p, x, OtherState, m, InitInp, ErrStat, ErrMsg )
   use YAML, only: yaml_write_var, yaml_write_array, yaml_write_list
      ! passed variables
   TYPE(BD_ParameterType),       INTENT(IN)     :: p                 !< Parameters of the structural dynamics module
   type(BD_ContinuousStateType), intent(in)     :: x                 !< Continuous states
   TYPE(BD_OtherStateType),      intent(in   ) :: OtherState         !< Other states at t
   TYPE(BD_MiscVarType),         INTENT(INout)  :: m                 !< misc/optimization variables
   TYPE(BD_InitInputType),       INTENT(IN   )  :: InitInp           !< Input data for initialization routine
   INTEGER(IntKi),               INTENT(OUT)    :: ErrStat           !< error status
   CHARACTER(*),                 INTENT(OUT)    :: ErrMsg            !< error message


      ! Local variables.

   INTEGER(IntKi)               :: I                                               ! Index for the nodes
   INTEGER(IntKi)               :: j, k                                            ! Generic index
   INTEGER(IntKi)               :: UnSu                                            ! I/O unit number for the summary output file

   CHARACTER(*), PARAMETER      :: FmtDat    = '(A,T35,1(:,F13.3))'                ! Format for outputting mass and modal data.
   CHARACTER(*), PARAMETER      :: FmtDatT   = '(A,T35,1(:,F13.8))'                ! Format for outputting time steps.

   CHARACTER(80)                :: OutPFmt                                         ! Format to print list of selected output channels to summary file

   ! Open the summary file and give it a heading.

   !$OMP critical(fileopen_critical)
   CALL GetNewUnit( UnSu, ErrStat, ErrMsg )
   CALL OpenFOutFile ( UnSu, TRIM( InitInp%RootName )//'.sum.yaml', ErrStat, ErrMsg )
   !$OMP end critical(fileopen_critical)
   IF ( ErrStat >= AbortErrLev ) RETURN

      ! Heading:
   WRITE (UnSu,'(A)')  '#This summary information was generated by '//TRIM( GetNVD(BeamDyn_Ver) )
                         
   WRITE (UnSu,'(/,A)')  '# --- Main parameters'
   call yaml_write_var  (UnSu, 'Mass'    , p%blade_mass      , 'F13.3', ErrStat, ErrMsg, comment='(kg)')
   call yaml_write_var  (UnSu, 'Length'  , p%blade_length    , 'F13.3', ErrStat, ErrMsg, comment='(m)')
   call yaml_write_list (UnSu, 'CG'      , p%blade_CG        , 'ES18.5', ErrStat, ErrMsg, comment='Blade center of mass (IEC coords) (m) from blade root')
   call yaml_write_array(UnSu, 'JRoot'   , p%blade_IN        , 'ES18.5', ErrStat, ErrMsg, comment='Blade mass moment of inertia at blade root. NOTE: from mass distribution only, missing some important inertial contributions (see PR#1337)')
   call yaml_write_list (UnSu, 'GlbPos'  , OtherState%GlbPos , 'ES18.5', ErrStat, ErrMsg, comment='Global position vector (IEC coords) of blade root at Initialization')
   call yaml_write_array(UnSu, 'GlbRot'  , OtherState%GlbRot , 'ES18.5', ErrStat, ErrMsg, comment='Global rotation tensor (IEC coords) at Initialization')
   call yaml_write_array(UnSu, 'RootOri' , InitInp%RootOri   , 'ES18.5', ErrStat, ErrMsg, comment='Initial blade orientation tensor (relative to global rotation tensor)')
   call yaml_write_list (UnSu, 'GlbCrv'  , OtherState%Glb_crv, 'ES18.5', ErrStat, ErrMsg, comment='Global rotation WM parameters (IEC coords) at Initialization')
   call yaml_write_list (UnSu, 'Gravity' , p%gravity         , 'ES18.5', ErrStat, ErrMsg, comment='Gravity vector (m/s^2) (IEC coords)')

!FIXME:analysis_type
   IF(p%analysis_type .EQ. BD_STATIC_ANALYSIS) THEN
       WRITE (UnSu,'(A,T59,A15)')  'Analysis_type:','"STATIC"'
   ELSEIF(p%analysis_type .EQ. BD_DYNAMIC_ANALYSIS) THEN
       WRITE (UnSu,'(A,T59,A15)')  'Analysis_type:','"DYNAMIC"'
   ELSEIF(p%analysis_type .EQ. BD_DYN_SSS_ANALYSIS) THEN
       WRITE (UnSu,'(A,T37,A)')  'Analysis_type:','"DYNAMIC with quasi-steady-state start"'
   ENDIF

   call yaml_write_var  (UnSu, 'Rhoinf'    , p%rhoinf    , 'ES15.5', ErrStat, ErrMsg, comment='Numerical damping parameter')
   call yaml_write_var  (UnSu, 'dt'        , p%dt        , 'ES15.5', ErrStat, ErrMsg, comment='Time increment')
   call yaml_write_var  (UnSu, 'niter'     , p%niter     , 'I15'    , ErrStat, ErrMsg, comment='Maximum number of iterations in Newton-Raphson solution')
   call yaml_write_var  (UnSu, 'tol'       , p%tol       , 'ES15.5' , ErrStat, ErrMsg, comment='Convergence parameter')
   call yaml_write_var  (UnSu, 'n_fact'    , p%n_fact    , 'I15'    , ErrStat, ErrMsg, comment='Factorization frequency in Newton-Raphson solution')

   IF(p%quadrature .EQ. GAUSS_QUADRATURE) THEN
       WRITE (UnSu,'(A,T59,A15)')  'Quadrature_method: ', 'Gaussian'
   ELSEIF(p%quadrature .EQ. TRAP_QUADRATURE) THEN
       WRITE (UnSu,'(A,T59,A15)')  'Quadrature_method: ', 'Trapezoidal'
       WRITE (UnSu,'(A,T59,I15)' ) 'FE_mesh_refinement_factor:', p%refine
   ENDIF
   WRITE (UnSu,'(A,T59,I15)' ) 'Number_of_elements:    ', p%elem_total
   WRITE (UnSu,'(A,T59,I15)' ) 'Number_of_nodes:       ', p%node_total

   WRITE (UnSu,'(/,A)')  '# --- Initial values and conditions'
   ! EB: NOTE: information about node and global node were lost here. Can be reintroduced by creating a matrix with indices and u in it. Or use "label" keyword in Yaml. 
   ! OK until we introduce more elements
   WRITE (UnSu,'(A)')  '# Initial position and Weiner-Milenkovic rotation vectors of nodes (IEC coordinate system)'
   WRITE (UnSu, '("#", 6(2x,A))') '            X        ','        Y        ','        Z        ','      WM_x       ','      WM_y       ','      WM_z       '
   WRITE (UnSu, '("#", 6(2x,A))') '    -----------------','-----------------','-----------------','-----------------','-----------------','-----------------'
   DO i=1,p%elem_total
       WRITE (UnSu, '("#", 1x,A,I4)') 'Element number: ',i
       call yaml_write_array(UnSu, 'Init_Nodes_E'//num2lstr(i), transpose(p%uuN0(1:6,:,i)), 'DummyFmt', ErrStat, ErrMsg, AllFmt='5(ES18.5,","),ES18.5')
   ENDDO

   WRITE (UnSu,'(/,A)')  '# Quadrature points position and rotation vectors'
   WRITE (UnSu, '("#", 6(2x,A))') '            X        ','        Y        ','        Z        ','      WM_x       ','      WM_y       ','      WM_z       '
   WRITE (UnSu, '("#", 6(2x,A))') '    -----------------','-----------------','-----------------','-----------------','-----------------','-----------------'
   DO i=1,p%elem_total
       WRITE (UnSu, '("#", 1x,A,I4)')  'Element number: ',i
       call yaml_write_array(UnSu, 'Init_QP_E'//num2lstr(i), transpose(p%uu0(1:6,:,i)), 'DummyFmt', ErrStat, ErrMsg, AllFmt='5(ES18.5,","),ES18.5')
   ENDDO

   WRITE (UnSu,'(/,A)')  '# Sectional stiffness and mass matrices at quadrature points (in IEC coordinates)'
   DO i=1,size(p%Stif0_QP,3)
      WRITE (UnSu,'(A,I4)')  '# Quadrature point number: ',i
      call yaml_write_array(UnSu, 'Init_K_QP'//num2lstr(i), p%Stif0_QP(:,1:6,i), 'ES15.5', ErrStat, ErrMsg)
      call yaml_write_array(UnSu, 'Init_M_QP'//num2lstr(i), p%Mass0_QP(:,1:6,i), 'ES15.5', ErrStat, ErrMsg)
   ENDDO
   call yaml_write_array(UnSu, 'Init_q'   , transpose(x%q(1:6,:)   ), 'ES18.5', ErrStat, ErrMsg, comment='Initial displacement and rotation')
   call yaml_write_array(UnSu, 'Init_dqdt', transpose(x%dqdt(1:6,:)), 'ES18.5', ErrStat, ErrMsg, comment='Initial velocity and angular velocity')


   WRITE (UnSu,'(/,A)')  '# --- Outputs'
   select case (p%BldMotionNodeLoc)
   case (BD_MESH_FE)  
      WRITE (UnSu, '(A)')  'Output_nodes_location: "finite element nodes"'
   case (BD_MESH_QP)  
      WRITE (UnSu, '(A)')  'Output_nodes_location: "quadrature points"'
   case (BD_MESH_STATIONS)  
      WRITE (UnSu, '(A)')  'Output_nodes_location: "blade propery input station locations"'
      !bjj: need to write where these nodes are located...
   end select
   

   
   ! Output channels:
   OutPFmt = '("# - [", I4, ",", 3X, A'//TRIM(Num2LStr(ChanLen+2))//', "," , 1 X, A'//TRIM(Num2LStr(ChanLen+2))//', "]" )'
   WRITE (UnSu, '(A)') "# OutList: # Requested Outputs (Index, Parameter, Unit)"
   DO I = 0,p%NumOuts
      WRITE (UnSu,OutPFmt)  I, '"'//trim(p%OutParam(I)%Name)//'"', '"'//trim(p%OutParam(I)%Units)//'"'
   END DO


   WRITE (UnSu,'("#", 15x,A)')
   WRITE (UnSu, '(A)') "# OutList_BldNd: # Requested Outputs at each blade station (Index, Parameter, Unit)"
   DO I = 1,p%BldNd_NumOuts
      WRITE (UnSu,OutPFmt)  I, '"'//trim(p%BldNd_OutParam(I)%Name)//'"', '"'//trim(p%BldNd_OutParam(I)%Units)//'"'
   END DO


   if ( p%analysis_type /= BD_STATIC_ANALYSIS ) then !dynamic analysis
      ! we'll add mass and stiffness matrices in the first call to UpdateStates
      m%Un_Sum = UnSu
      WRITE (UnSu,'(/,A)')  '# --- System matrices'
   else
      CLOSE(UnSu)
   end if


RETURN
END SUBROUTINE BD_PrintSum
!----------------------------------------------------------------------------------------------------------------------------------

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!> This routine initializes the array that maps rows/columns of the Jacobian to specific mesh fields.
!! Do not change the order of this packing without changing subroutine !
SUBROUTINE Init_Jacobian( p, u, y, m, InitOut, ErrStat, ErrMsg)

   TYPE(BD_ParameterType)            , INTENT(INOUT) :: p                     !< parameters
   TYPE(BD_InputType)                , INTENT(IN   ) :: u                     !< inputs
   TYPE(BD_OutputType)               , INTENT(IN   ) :: y                     !< outputs
   TYPE(BD_MiscVarType)              , INTENT(INOUT) :: m                     !< misc var data
   TYPE(BD_InitOutputType)           , INTENT(INOUT) :: InitOut               !< Initialization output data (for Jacobian row/column names)
   
   INTEGER(IntKi)                    , INTENT(  OUT) :: ErrStat               !< Error status of the operation
   CHARACTER(*)                      , INTENT(  OUT) :: ErrMsg                !< Error message if ErrStat /= ErrID_None
   
   INTEGER(IntKi)                                    :: ErrStat2
   CHARACTER(ErrMsgLen)                              :: ErrMsg2
   CHARACTER(*), PARAMETER                           :: RoutineName = 'Init_Jacobian'
   
      ! local variables:
   INTEGER(IntKi)                :: i, j, index, nu, i_meshField
   REAL(R8Ki)                    :: perturb, perturb_b
   REAL(R8Ki)                    :: MaxThrust, MaxTorque
   CHARACTER(1), PARAMETER       :: UVW(3) = (/'U','V','W'/)
   
            
   
   ErrStat = ErrID_None
   ErrMsg  = ""
   
   call Init_Jacobian_y( p, y, InitOut, ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)

   call Init_Jacobian_x_z( p, InitOut, ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   if (ErrStat >= AbortErrLev) return

      ! determine how many inputs there are in the Jacobians
   if (p%CompAeroMaps) then
      nu = u%DistrLoad%NNodes  * 6               ! 3 forces + 3 moments at each node
   else
      nu = u%RootMotion%NNodes * 18            & ! 3 Translation Displacements + 3 orientations + 6 velocities (rotation+translation) + 6 accelerations at each node
         + u%PointLoad%NNodes  * 6             & ! 3 forces + 3 moments at each node
         + u%DistrLoad%NNodes  * 6               ! 3 forces + 3 moments at each node
   end if
   
   ! all other inputs (e.g., hub motion) ignored

   !............................
   ! fill matrix to store index to help us figure out what the ith value of the u vector really means
   ! (see beamdyn::perturb_u ... these MUST match )
   ! column 1 indicates module's mesh and field
   ! column 2 indicates the first index (x-y-z component) of the field
   ! column 3 is the node
   !............................
   
   call allocAry( p%Jac_u_indx, nu, 3, 'p%Jac_u_indx', ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   if (ErrStat >= AbortErrLev) return
            
   !...............
   ! BD input mappings stored in p%Jac_u_indx:   
   !...............            
   index = 1
   !Module/Mesh/Field: u%RootMotion%TranslationDisp  = 1;
   !Module/Mesh/Field: u%RootMotion%Orientation      = 2;
   !Module/Mesh/Field: u%RootMotion%TranslationVel   = 3;
   !Module/Mesh/Field: u%RootMotion%RotationVel      = 4;
   !Module/Mesh/Field: u%RootMotion%TranslationAcc   = 5;
   !Module/Mesh/Field: u%RootMotion%RotationAcc      = 6;
   if (.not. p%CompAeroMaps) then
      do i_meshField = 1,6
         do i=1,u%RootMotion%NNodes
            do j=1,3
               p%Jac_u_indx(index,1) =  i_meshField
               p%Jac_u_indx(index,2) =  j !component index:  j
               p%Jac_u_indx(index,3) =  i !Node:   i
               index = index + 1
            end do !j
         end do !i
      end do
   
      !Module/Mesh/Field: u%PointLoad%Force   = 7;
      !Module/Mesh/Field: u%PointLoad%Moment  = 8;
      do i_meshField = 7,8
         do i=1,u%PointLoad%NNodes
            do j=1,3
               p%Jac_u_indx(index,1) =  i_meshField
               p%Jac_u_indx(index,2) =  j !component index:  j
               p%Jac_u_indx(index,3) =  i !Node:   i
               index = index + 1
            end do !j
         end do !i
      end do
   end if
   
   !Module/Mesh/Field: u%DistrLoad%Force   =  9;
   !Module/Mesh/Field: u%DistrLoad%Moment  = 10;
   do i_meshField = 9,10
      do i=1,u%DistrLoad%NNodes
         do j=1,3
            p%Jac_u_indx(index,1) =  i_meshField
            p%Jac_u_indx(index,2) =  j !component index:  j
            p%Jac_u_indx(index,3) =  i !Node:   i
            index = index + 1
         end do !j
      end do !i
   end do !i_meshField
      
   
   
      !......................................
      ! default perturbations, p%du:
      !......................................
   call allocAry( p%du, 10, 'p%du', ErrStat2, ErrMsg2) ! 10 = number of unique values in p%Jac_u_indx(:,1)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)

   perturb   = 0.2_R8Ki*D2R_D
   perturb_b = 0.2_R8Ki*D2R_D * p%blade_length
   
   MaxThrust = 170.0_R8Ki*p%blade_length**2
   MaxTorque =  14.0_R8Ki*p%blade_length**3
   
   p%du( 1) = perturb_b                                        ! u%RootMotion%TranslationDisp  = 1;
   p%du( 2) = perturb                                          ! u%RootMotion%Orientation      = 2;
   p%du( 3) = perturb_b                                        ! u%RootMotion%TranslationVel   = 3;
   p%du( 4) = perturb                                          ! u%RootMotion%RotationVel      = 4;
   p%du( 5) = perturb_b                                        ! u%RootMotion%TranslationAcc   = 5;
   p%du( 6) = perturb                                          ! u%RootMotion%RotationAcc      = 6;

   p%du( 7) = MaxThrust / (100.0_R8Ki * 3.0_R8Ki * u%PointLoad%NNodes )  ! u%PointLoad%Force             = 7;
   p%du( 8) = MaxTorque / (100.0_R8Ki * 3.0_R8Ki * u%PointLoad%NNodes )  ! u%PointLoad%Moment            = 8;
   
   p%du( 9) = MaxThrust / (100.0_R8Ki * 3.0_R8Ki * u%DistrLoad%NNodes )  ! u%DistrLoad%Force             = 9;
   p%du(10) = MaxTorque / (100.0_R8Ki * 3.0_R8Ki * u%DistrLoad%NNodes )  ! u%DistrLoad%Moment            =10;
   
      !.....................
      ! get names of linearized inputs
      !.....................
   call AllocAry(InitOut%LinNames_u, nu, 'LinNames_u', ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   call AllocAry(InitOut%RotFrame_u, nu, 'RotFrame_u', ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   call AllocAry(InitOut%IsLoad_u, nu, 'IsLoad_u', ErrStat2, ErrMsg2)
      call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
      if (ErrStat >= AbortErrLev) return

   InitOut%RotFrame_u = .false. ! every input is on a mesh, which stores values in the global (not rotating) frame
   
   index = 1
   InitOut%IsLoad_u = .true.  ! initialize all inputs as loads, and overwrite for the RootMotion mesh, below:
   if (.not. p%CompAeroMaps) then
      call PackMotionMesh_Names(u%RootMotion, 'RootMotion', InitOut%LinNames_u, index) ! all 6 motion fields
      InitOut%IsLoad_u(1:index-1) = .false. ! the RootMotion inputs are not loads
      call PackLoadMesh_Names(  u%PointLoad,   'PointLoad', InitOut%LinNames_u, index)
   end if
   call PackLoadMesh_Names(  u%DistrLoad,   'DistrLoad', InitOut%LinNames_u, index)


END SUBROUTINE Init_Jacobian
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine initializes the Jacobian parameters and initialization outputs for the linearized outputs.
SUBROUTINE Init_Jacobian_y( p, y, InitOut, ErrStat, ErrMsg)

   TYPE(BD_ParameterType)            , INTENT(INOUT) :: p                     !< parameters
   TYPE(BD_OutputType)               , INTENT(IN   ) :: y                     !< outputs
   TYPE(BD_InitOutputType)           , INTENT(INOUT) :: InitOut               !< Initialization output data (for Jacobian row/column names)
   
   INTEGER(IntKi)                    , INTENT(  OUT) :: ErrStat               !< Error status of the operation
   CHARACTER(*)                      , INTENT(  OUT) :: ErrMsg                !< Error message if ErrStat /= ErrID_None
   
      ! local variables:
   INTEGER(IntKi)                                    :: i
   INTEGER(IntKi)                                    :: j
   INTEGER(IntKi)                                    :: index_next
   LOGICAL                                           :: AllOut(MaxOutPts)
   INTEGER(IntKi)                                    :: ErrStat2
   CHARACTER(ErrMsgLen)                              :: ErrMsg2
   CHARACTER(*), PARAMETER                           :: RoutineName = 'Init_Jacobian_y'

   CHARACTER(ChanLen)                                :: ChannelName
   LOGICAL                                           :: isRotating
   LOGICAL                                           :: BladeMask(FIELDMASK_SIZE)   ! flags to determine if this field is part of the packing
   
   ErrStat = ErrID_None
   ErrMsg  = ""
   
   if (p%CompAeroMaps) then
      p%Jac_ny = y%BldMotion%NNodes     * 12     ! 6 displacements (translation, rotation) + 6 velocities 
   else
   
         ! determine how many outputs there are in the Jacobians     
      p%Jac_ny = y%ReactionForce%NNodes * 6     & ! 3 forces + 3 moments at each node
               + y%BldMotion%NNodes     * 18    & ! 6 displacements (translation, rotation) + 6 velocities + 6 accelerations at each node
               + p%NumOuts + p%BldNd_TotNumOuts   ! WriteOutput values 
   end if
   
      ! get the names of the linearized outputs:
   call AllocAry(InitOut%LinNames_y, p%Jac_ny,'LinNames_y',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
   call AllocAry(InitOut%RotFrame_y, p%Jac_ny,'RotFrame_y',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
      if (ErrStat >= AbortErrLev) return
   
         
   InitOut%RotFrame_y = .false.   ! need to set all the values in the global system to .false
   
   index_next = 1
   if (p%CompAeroMaps) then
      BladeMask = .true. ! default is all the fields
      BladeMask(MASKID_TRANSLATIONACC) = .false.
      BladeMask(MASKID_ROTATIONACC)    = .false.
   
      call PackMotionMesh_Names(y%BldMotion,     'Blade motion',   InitOut%LinNames_y, index_next, FieldMask=BladeMask)
   else
      call PackLoadMesh_Names(  y%ReactionForce, 'Reaction force', InitOut%LinNames_y, index_next)
      call PackMotionMesh_Names(y%BldMotion,     'Blade motion',   InitOut%LinNames_y, index_next)

      do i=1,p%NumOuts + p%BldNd_TotNumOuts
         InitOut%LinNames_y(i+index_next-1) = trim(InitOut%WriteOutputHdr(i))//', '//trim(InitOut%WriteOutputUnt(i))
      end do
   
      AllOut = .true. ! all output values except those specifically in the global system are in the rotating system
      AllOut(TipTVXg) = .false.
      AllOut(TipTVYg) = .false.
      AllOut(TipTVZg) = .false.
      AllOut(TipRVXg) = .false.
      AllOut(TipRVYg) = .false.
      AllOut(TipRVZg) = .false.
      
      do j=1,9
         do i=1,3 !x,y,z
            AllOut(NTVg(j,i)) = .false.
            AllOut(NRVg(j,i)) = .false.
         end do
      end do
   
      do i=1,p%NumOuts
         if (p%OutParam(i)%Indx == 0 ) then
            InitOut%RotFrame_y(i+index_next-1) = .false.
         else
            InitOut%RotFrame_y(i+index_next-1) = AllOut( p%OutParam(i)%Indx )
         end if
      end do
   

         ! set outputs for all nodes out:
      index_next = index_next + p%NumOuts
      DO i=1,p%BldNd_NumOuts
         ChannelName = p%BldNd_OutParam(i)%Name
         call Conv2UC(ChannelName)
         if ( ChannelName( LEN_TRIM(ChannelName):LEN_TRIM(ChannelName) ) == 'G') then ! channel is in global coordinate system
            isRotating = .false.
         else
            isRotating = .true.
         end if
         InitOut%RotFrame_y(index_next : index_next+size(p%BldNd_BlOutNd)-1 ) = isRotating
         index_next = index_next + size(p%BldNd_BlOutNd)
      ENDDO
   end if


END SUBROUTINE Init_Jacobian_y
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine initializes the Jacobian parameters and initialization outputs for the linearized continuous states.
SUBROUTINE Init_Jacobian_x_z( p, InitOut, ErrStat, ErrMsg)

   TYPE(BD_ParameterType)            , INTENT(INOUT) :: p                     !< parameters
   TYPE(BD_InitOutputType)           , INTENT(INOUT) :: InitOut               !< Output for initialization routine   
   
   INTEGER(IntKi)                    , INTENT(  OUT) :: ErrStat               !< Error status of the operation
   CHARACTER(*)                      , INTENT(  OUT) :: ErrMsg                !< Error message if ErrStat /= ErrID_None
   
   INTEGER(IntKi)                                    :: ErrStat2
   CHARACTER(ErrMsgLen)                              :: ErrMsg2
   CHARACTER(*), PARAMETER                           :: RoutineName = 'Init_Jacobian_x'
   CHARACTER(200)                                    :: Describe
   
      ! local variables:
   INTEGER(IntKi)                :: i
   INTEGER(IntKi)                :: indx
   
   ErrStat = ErrID_None
   ErrMsg  = ""
   
   p%Jac_nx = p%dof_node * (p%node_total-1) ! the first node is actually a constraint state

      ! allocate space for the row/column names and for perturbation sizes
   !call allocAry(p%dx,               p%dof_node*(p%node_total-1),     'p%dx',       ErrStat2, ErrMsg2); call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   CALL AllocAry(InitOut%LinNames_x, p%Jac_nx*2, 'LinNames_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   CALL AllocAry(InitOut%RotFrame_x, p%Jac_nx*2, 'RotFrame_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   CALL AllocAry(InitOut%DerivOrder_x, p%Jac_nx*2, 'DerivOrder_x', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   !CALL AllocAry(InitOut%LinNames_z, p%dof_node*2, 'LinNames_z', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   !CALL AllocAry(InitOut%RotFrame_z, p%dof_node*2, 'RotFrame_z', ErrStat2, ErrMsg2); CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
   if (ErrStat >= AbortErrLev) return
   
   
      !......................................
      ! default perturbations, p%dx:
      !......................................
   p%dx(1:3) = 0.2_BDKi*D2R_D * p%blade_length                    ! deflection states in m and m/s
   p%dx(4:6) = 0.2_BDKi*D2R_D                                     ! deflection states in rad and rad/s
   
   InitOut%RotFrame_x   = p%RotStates
   InitOut%DerivOrder_x = 2
   
      !......................................
      ! set linearization output names:
      !......................................
   indx = 1
   DO i=2, p%node_total
      Describe = 'finite element node '//trim(num2lstr(i))//' (number of elements = '//trim(num2lstr(p%elem_total))//'; element order = '//trim(num2lstr(p%nodes_per_elem-1))//')'
      InitOut%LinNames_x(indx) = trim(Describe)//' translational displacement in X, m'
      indx = indx + 1
      InitOut%LinNames_x(indx) = trim(Describe)//' translational displacement in Y, m'
      indx = indx + 1
      InitOut%LinNames_x(indx) = trim(Describe)//' translational displacement in Z, m'
      indx = indx + 1
      InitOut%LinNames_x(indx) = trim(Describe)//' rotational displacement in X, rad'
      indx = indx + 1
      InitOut%LinNames_x(indx) = trim(Describe)//' rotational displacement in Y, rad'
      indx = indx + 1
      InitOut%LinNames_x(indx) = trim(Describe)//' rotational displacement in Z, rad'
      indx = indx + 1
   END DO
   
   do i=1,p%Jac_nx
      InitOut%LinNames_x(i+p%Jac_nx) = 'First time derivative of '//trim(InitOut%LinNames_x(i))//'/s'
      InitOut%RotFrame_x(i+p%Jac_nx) = InitOut%RotFrame_x(i)
   end do


   !InitOut%RotFrame_z    = .true.
   !InitOut%LinNames_z(1) = 'Node 1 translational displacement in X, m'
   !InitOut%LinNames_z(2) = 'Node 1 translational displacement in Y, m'
   !InitOut%LinNames_z(3) = 'Node 1 translational displacement in Z, m'
   !InitOut%LinNames_z(4) = 'Node 1 rotational displacement in X, -'
   !InitOut%LinNames_z(5) = 'Node 1 rotational displacement in Y, -'
   !InitOut%LinNames_z(6) = 'Node 1 rotational displacement in Z, -'
   !
   !do i=1,6
   !   InitOut%LinNames_x(i+6) = 'First time derivative of '//trim(InitOut%LinNames_z(i))//'/s'
   !end do
   
   
END SUBROUTINE Init_Jacobian_x_z
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine perturbs the nth element of the u array (and mesh/field it corresponds to)
!! Do not change this without making sure subroutine beamdyn::init_jacobian is consistant with this routine!
SUBROUTINE Perturb_u( p, n, perturb_sign, u, du )

   TYPE(BD_ParameterType)              , INTENT(IN   ) :: p                      !< parameters
   INTEGER( IntKi )                    , INTENT(IN   ) :: n                      !< number of array element to use 
   INTEGER( IntKi )                    , INTENT(IN   ) :: perturb_sign           !< +1 or -1 (value to multiply perturbation by; positive or negative difference)
   TYPE(BD_InputType)                  , INTENT(INOUT) :: u                      !< perturbed BD inputs
   REAL( R8Ki )                        , INTENT(  OUT) :: du                     !< amount that specific input was perturbed
   

   ! local variables
   INTEGER                                             :: fieldIndx
   INTEGER                                             :: node
      
   fieldIndx = p%Jac_u_indx(n,2)
   node      = p%Jac_u_indx(n,3)
   
   du = p%du(  p%Jac_u_indx(n,1) )
   
      ! determine which mesh we're trying to perturb and perturb the input:
   SELECT CASE( p%Jac_u_indx(n,1) )
      
   CASE ( 1) !Module/Mesh/Field: u%RootMotion%TranslationDisp = 1;
      u%RootMotion%TranslationDisp( fieldIndx,node) = u%RootMotion%TranslationDisp( fieldIndx,node) + du * perturb_sign
   CASE ( 2) !Module/Mesh/Field: u%RootMotion%Orientation = 2;
      CALL PerturbOrientationMatrix( u%RootMotion%Orientation(:,:,node), du * perturb_sign, fieldIndx )   ! NOTE: call not using DCM_logmap
   CASE ( 3) !Module/Mesh/Field: u%RootMotion%TranslationVel = 3;
      u%RootMotion%TranslationVel( fieldIndx,node) = u%RootMotion%TranslationVel( fieldIndx,node) + du * perturb_sign
   CASE ( 4) !Module/Mesh/Field: u%RootMotion%RotationVel = 4;
      u%RootMotion%RotationVel(fieldIndx,node) = u%RootMotion%RotationVel(fieldIndx,node) + du * perturb_sign
   CASE ( 5) !Module/Mesh/Field: u%RootMotion%TranslationAcc = 5;
      u%RootMotion%TranslationAcc( fieldIndx,node) = u%RootMotion%TranslationAcc( fieldIndx,node) + du * perturb_sign
   CASE ( 6) !Module/Mesh/Field: u%RootMotion%RotationAcc = 6;
      u%RootMotion%RotationAcc(fieldIndx,node) = u%RootMotion%RotationAcc(fieldIndx,node) + du * perturb_sign
   
   CASE ( 7) !Module/Mesh/Field: u%PointLoad%Force = 7;
      u%PointLoad%Force(fieldIndx,node) = u%PointLoad%Force(fieldIndx,node) + du * perturb_sign 
   CASE ( 8) !Module/Mesh/Field: u%PointLoad%Moment  = 8;
      u%PointLoad%Moment(fieldIndx,node) = u%PointLoad%Moment(fieldIndx,node) + du * perturb_sign
      
   CASE ( 9) !Module/Mesh/Field: u%DistrLoad%Force   =  9;
      u%DistrLoad%Force( fieldIndx,node) = u%DistrLoad%Force( fieldIndx,node) + du * perturb_sign
   CASE (10) !Module/Mesh/Field: u%DistrLoad%Moment = 10;
      u%DistrLoad%Moment(fieldIndx,node) = u%DistrLoad%Moment(fieldIndx,node) + du * perturb_sign
      
   END SELECT
      
END SUBROUTINE Perturb_u
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine uses values of two output types to compute an array of differences.
!! Do not change this packing without making sure subroutine beamdyn::init_jacobian is consistant with this routine!
SUBROUTINE Compute_dY(p, y_p, y_m, delta, dY)
   
   TYPE(BD_ParameterType)            , INTENT(IN   ) :: p         !< parameters
   TYPE(BD_OutputType)               , INTENT(IN   ) :: y_p       !< BD outputs at \f$ u + \Delta_p u \f$ or \f$ z + \Delta_p z \f$ (p=plus)
   TYPE(BD_OutputType)               , INTENT(IN   ) :: y_m       !< BD outputs at \f$ u - \Delta_m u \f$ or \f$ z - \Delta_m z \f$ (m=minus)   
   REAL(R8Ki)                        , INTENT(IN   ) :: delta     !< difference in inputs or states \f$ delta_p = \Delta_p u \f$ or \f$ delta_p = \Delta_p x \f$
   REAL(R8Ki)                        , INTENT(INOUT) :: dY(:)     !< column of dYdu or dYdx: \f$ \frac{\partial Y}{\partial u_i} = \frac{y_p - y_m}{2 \, \Delta u}\f$ or \f$ \frac{\partial Y}{\partial z_i} = \frac{y_p - y_m}{2 \, \Delta x}\f$
   
      ! local variables:
   INTEGER(IntKi)                                    :: i              ! loop over outputs
   INTEGER(IntKi)                                    :: indx_first     ! index indicating next value of dY to be filled 
   LOGICAL                                           :: Mask(FIELDMASK_SIZE)   ! flags to determine if this field is part of the packing

   indx_first = 1
   if (p%CompAeroMaps) then
      Mask  = .true.
      Mask(MASKID_TRANSLATIONACC) = .false.
      Mask(MASKID_ROTATIONACC) = .false.
      call PackMotionMesh_dY(y_p%BldMotion,     y_m%BldMotion,     dY, indx_first, FieldMask=Mask) ! 4 motion fields
   else
      call PackLoadMesh_dY(  y_p%ReactionForce, y_m%ReactionForce, dY, indx_first)
      call PackMotionMesh_dY(y_p%BldMotion,     y_m%BldMotion,     dY, indx_first) ! all 6 motion fields
   
      do i=1,p%NumOuts + p%BldNd_TotNumOuts
         dY(i+indx_first-1) = y_p%WriteOutput(i) - y_m%WriteOutput(i)
      end do
   end if
   
   
   dY = dY / (2.0_R8Ki*delta)
   
END SUBROUTINE Compute_dY
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine perturbs the nth element of the x array (and mesh/field it corresponds to)
!! Do not change this without making sure subroutine beamdyn::init_jacobian is consistant with this routine!
SUBROUTINE Perturb_x( p, fieldIndx, node, dof, perturb_sign, x, dx )

   TYPE(BD_ParameterType)              , INTENT(IN   ) :: p                      !< parameters
   INTEGER( IntKi )                    , INTENT(IN   ) :: fieldIndx              !< field in the state type: 1=displacements; 2=velocities
   INTEGER( IntKi )                    , INTENT(IN   ) :: node                   !< node number
   INTEGER( IntKi )                    , INTENT(IN   ) :: dof                    !< dof for this perturbation
   INTEGER( IntKi )                    , INTENT(IN   ) :: perturb_sign           !< +1 or -1 (value to multiply perturbation by; positive or negative difference)
   TYPE(BD_ContinuousStateType)        , INTENT(INOUT) :: x                      !< perturbed BD states
   REAL( R8Ki )                        , INTENT(  OUT) :: dx                     !< amount that specific state was perturbed
   

   ! local variables
   integer(intKi)                                      :: ErrStat2
   character(ErrMsgLen)                                :: ErrMsg2
   
   REAL(R8Ki)                                          :: orientation(3,3)
   REAL(R8Ki)                                          :: rotation(3,3)
   REAL(R8Ki)                                          :: CrvPerturb(3), CrvBase(3)
   
   dx = p%dx(dof)
               
   if (fieldIndx==1) then
      if (dof < 4) then ! translational displacement
         x%q( dof, node ) = x%q( dof, node ) + dx * perturb_sign
      else ! w-m parameters
         
         ! Calculate perturbation in WM parameters
         CrvPerturb = 0.0_R8Ki
         CrvPerturb(dof-3) = 4.0_R8Ki * tan(dx * perturb_sign / 4.0_R8Ki)

         ! Get base rotation in WM parameters
         CrvBase = x%q(4:6, node)

         ! Compose pertubation and base rotation and store in state
         call BD_CrvCompose(x%q(4:6, node), CrvPerturb, CrvBase, FLAG_R1R2)
      end if
   else
      x%dqdt( dof, node ) = x%dqdt( dof, node ) + dx * perturb_sign
   end if

      
END SUBROUTINE Perturb_x
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine uses values of two output types to compute an array of differences.
!! Do not change this packing without making sure subroutine beamdyn::init_jacobian is consistant with this routine!
SUBROUTINE Compute_dX(p, x_p, x_m, delta, dX)
   
   TYPE(BD_ParameterType)            , INTENT(IN   ) :: p         !< parameters
   TYPE(BD_ContinuousStateType)      , INTENT(IN   ) :: x_p       !< BD continuous states at \f$ u + \Delta_p u \f$ or \f$ x + \Delta_p x \f$ (p=plus)
   TYPE(BD_ContinuousStateType)      , INTENT(IN   ) :: x_m       !< BD continuous states at \f$ u - \Delta_m u \f$ or \f$ x - \Delta_m x \f$ (m=minus)   
   REAL(R8Ki)                        , INTENT(IN   ) :: delta     !< difference in inputs or states \f$ delta_p = \Delta_p u \f$ or \f$ delta_p = \Delta_p x \f$
   REAL(R8Ki)                        , INTENT(INOUT) :: dX(:)     !< column of dXdu or dXdx: \f$ \frac{\partial X}{\partial u_i} = \frac{x_p - x_m}{2 \, \Delta u}\f$ or \f$ \frac{\partial X}{\partial x_i} = \frac{x_p - x_m}{2 \, \Delta x}\f$
   
      ! local variables:
   INTEGER(IntKi)                                    :: i              ! loop over nodes
   INTEGER(IntKi)                                    :: dof            ! loop over dofs
   INTEGER(IntKi)                                    :: index          ! index indicating next value of dX to be filled 
   
   index = 1
   do i=2,p%node_total
      do dof=1,p%dof_node
         dX(index) = x_p%q( dof, i ) - x_m%q( dof, i )
         index = index+1
      end do
   end do

   do i=2,p%node_total
      do dof=1,p%dof_node
         dX(index) = x_p%dqdt( dof, i ) - x_m%dqdt( dof, i )
         index = index+1
      end do
   end do

   dX = dX / ( 2.0_R8Ki*delta)
   
END SUBROUTINE Compute_dX
!----------------------------------------------------------------------------------------------------------------------------------
!> This routine uses values of two output types to compute an array of differences.
!! Do not change this packing without making sure subroutine beamdyn::init_jacobian is consistant with this routine!
SUBROUTINE Compute_RelState_Matrix(p, u, x, OtherState, RelState_x, RelState_xdot)
   
   TYPE(BD_ParameterType)            , INTENT(IN   ) :: p                  !< parameters
   TYPE(BD_InputType)                , INTENT(IN   ) :: u                  !< BD inputs
   TYPE(BD_ContinuousStateType)      , INTENT(IN   ) :: x                  !< BD continuous states
   TYPE(BD_OtherStateType)           , INTENT(IN   ) :: OtherState         !< Other states at t
   REAL(R8Ki)                        , INTENT(INOUT) :: RelState_x(:,:)    !< 
   REAL(R8Ki)                        , INTENT(INOUT) :: RelState_xdot(:,:) !< 
   
      ! local variables:
   INTEGER(IntKi)                                    :: i              ! loop counter
   INTEGER(IntKi)                                    :: j              ! loop counter
   INTEGER(IntKi)                                    :: dof            ! loop over dofs
   INTEGER(IntKi)                                    :: q_index        ! index into the state arrays
   INTEGER(IntKi)                                    :: dqdt_index     ! index into the state arrays
   INTEGER(IntKi)                                    :: node           ! node in the state arrays
   
   REAL(R8Ki)                                        :: dp             ! temporary dot product
   REAL(R8Ki)                                        :: cp(3)          ! temporary cross product
   REAL(R8Ki)                                        :: RotVel(3)      ! temporary velocity
   REAL(R8Ki)                                        :: RotAcc(3)      ! temporary acceleration
   REAL(R8Ki)                                        :: DisplacedPosition(3)
   REAL(R8Ki)                                        :: fx_p(3,3)
   
   RelState_x    = 0.0_ReKi
   RelState_xdot = 0.0_ReKi
   
   !-----------------------------------
   do i=1,p%elem_total
      do j=2,p%nodes_per_elem

         node       = (i-1)*(p%nodes_per_elem-1) + j  ! index to state array (rows of conversion matrices)
         q_index    = (node - 2)*p%dof_node + 1       ! index into displacement portion of x (skipping node 1)
         dqdt_index = p%Jac_nx + q_index
         
         DisplacedPosition = u%RootMotion%Position(:,1) + u%RootMotion%TranslationDisp(:,1) &
                           - OtherState%GlbPos - MATMUL(OtherState%GlbRot, p%uuN0(1:3,j,i) + x%q(1:3,node) )

         RotVel = real(u%RootMotion%RotationVel(:,1),R8Ki)
         RotAcc = real(u%RootMotion%RotationAcc(:,1),R8Ki)

         fx_p = SkewSymMat(DisplacedPosition)
         
         do dof=0,5
            RelState_x( q_index+dof, 1+dof   ) = 1.0_R8Ki      ! root displacements to node displacements
         end do
         do dof=0,5
            RelState_x( dqdt_index+dof, 7+dof   ) = 1.0_R8Ki   ! root velocities to node velocities
         end do
         
         
         RelState_x( q_index:q_index+2,        4: 6 ) = fx_p                                                            ! root rotational displacement to node translational displacement
         RelState_x( dqdt_index:dqdt_index+2, 10:12 ) = fx_p                                                            ! root rotational velocity to node translational velocity

            ! root rotational displacement to node translational velocity:
         RelState_x( dqdt_index:dqdt_index+2, 4:6   ) = OuterProduct( DisplacedPosition, RotVel )
         dp = dot_product( DisplacedPosition, RotVel )
         do dof=0,2
            RelState_x( dqdt_index+dof, 4+dof   ) = RelState_x( dqdt_index+dof, 4+dof   ) - dp                          ! root rotational displacement to node translational velocity 
         end do
         !----------
         
         
         !.............................................
         ! The first p%Jac_nx rows of RelState_xdot are the same as the last p%Jac_nx rows of RelState_x, so I'm not going to recalculate these rows, we'll set them after the loops:
         !do dof=0,5
         !   RelState_xdot( q_index+dof, 7+dof   ) = 1.0_ReKi      ! root velocities to node velocities
         !end do
         !RelState_xdot( q_index:q_index+2, 4:6 )         = RelState_x( dqdt_index:dqdt_index+2, 4:6 )                  ! root rotational displacement to node translational velocity
         !RelState_xdot( q_index:q_index+2,       10:12 ) = fx_p                                                        ! root rotational velocity to node translational velocity
         
         do dof=0,5
            RelState_xdot( dqdt_index+dof, 13+dof   ) = 1.0_R8Ki   ! root accelerations to node accelerations
         end do
         
         
            ! root translational velocity to node translational acceleration:
         cp = cross_product(u%RootMotion%RotationVel(:,1), DisplacedPosition)
         RelState_xdot( dqdt_index:dqdt_index+2, 7:9   ) = OuterProduct( DisplacedPosition, RotAcc )  &
                                                         + OuterProduct( cp, RotVel ) - dp*SkewSymMat(RotVel)
         dp = dot_product( DisplacedPosition, RotAcc )
         do dof=0,2
            RelState_xdot( dqdt_index+dof, 7+dof   ) = RelState_xdot( dqdt_index+dof, 7+dof   ) - dp
         end do
         !-----------
         
         RelState_xdot( dqdt_index:dqdt_index+2, 10:12 ) = RelState_x( dqdt_index:dqdt_index+2, 4:6 ) + SkewSymMat(cp)  ! root rotational velocity to node translational acceleration
         RelState_xdot( dqdt_index:dqdt_index+2, 16:18 ) = fx_p                                                         ! root rotational acceleration to node translational acceleration
         
      end do
   end do
   RelState_xdot(1:p%Jac_nx,:) = RelState_x(p%Jac_nx+1:,:)

END SUBROUTINE Compute_RelState_Matrix
!----------------------------------------------------------------------------------------------------------------------------------

!----------------------------------------------------------------------------------------------------------------------------------
END MODULE BeamDyn_IO
